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CHAPTER  I 


INTRODUCTION 


Discrete  data  sampling  at  the  Nyquist  rate  (or  better)  for  the 
reconstruction  of  continuous  waveforms  is  well  known.  However,  the 
application  of  the  sampling  theorem  [1,2]  is  often  limited  to  one 
dimensional  problems.  This  report  will  explore  an  application  of  the  N 
dimensional  sampling  theorem  to  inverse  scattering.  More  specifically, 
the  application  is  on  the  reconstruction  of  the  spatial  impulse  response 
of  a  finite  object  using  different  sampling  criteria.  The  different 
aspects  of  the  N  dimensional  sampling  theorem  are  investigated  to 
produce  efficient  sampling  criteria  in  the  wave  number  space  such  that 
the  spatial  impulse  response  of  a  finite  object  is  sufficiently 
characterized.  In  addition,  the  impulse  response  concept  is  extended  to 
create  possibly  an  image  of  the  object. 

The  impulse  response  [3,4]  is  important  both  in  target 
identification  and  imaging.  The  impulse  response  concept,  as  most 
people  understand,  is  a  far  field  one  dimensional  time  response  concept. 
This  time  response  waveform  when  applied  in  scattering  can  be  plotted  on 
a  distance  abscissa  after  the  factor  due  to  the  speed  of  the  wave  is 
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n 
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accounted  for.  If  one  can  obtain  the  frequency  responses  of  a  finite 
object  for  all  aspect  angles  over  the  4*  solid  angle,  the  inverse 
Fourier  transform  of  the  total  frequency  response  into  the  spatial 
domain  forms  an  image  which  is  called  the  spatial  impulse  response, 
i.e.,  the  impulse  response  of  a  finite  object  at  x  is 


A  1  °° 

f<;>  “727)7  / 


-  j (C-*)  . 

F(k)e  dk 


(1-1) 


where  F(k)  =  the  far  field  frequency  response  of  a  finite 
object  at  k  in  the  space-frequency  domain 

x  "  ^xl»  x2*  x3^ 
k  =  (kp  k2,  k3) 

dk  =  dkx  dk2  dk3 

The  three  dimensional  function  f(x)  for  all  x's  in  the  x-space 
constitutes  an  image  called  the  spatial  impulse  response. 

The  spatial  impluse  response  can  be  divided  into  two  types  -  the 
monostatic  and  the  bistatic.  The  monostatic  spatial  impulse  response 
uses  monostatic  frequency  responses  of  all  aspect  angles.  The  bistatic 
spatial  impulse  response  uses  the  bistatic  frequency  responses  of  all 
receiving  aspect  angles.  The  spatial,  two  dimensional,  or  one 
dimensional  impulse  responses  discussed  hereafter  wi 1 1  all  be  referring 
to  the  monostatic  case  unless  otherwise  stated.  The  data  used  are  also 
monostatic  data.  However,  the  theory  to  be  discussed  is  also  applicable 
to  bistatic  cases.  Furthermore,  only  two  dimensional  frequency  data  are 
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used  in  the  discussion,  so  this  report  will  emphasize  the  two 
dimensional  impulse  response  only. 

In  the  introduction,  the  theme  and  purpose  of  this  report  are 
defined.  A  brief  description  on  each  of  the  following  chapters  is  also 
included.  In  Chapter  II,  the  impulse  response,  sampling  theorem  and 
their  relationship  are  discussed.  The  basics  of  the  one  dimensional 
impulse  response  concept  and  the  one  dimensional  sampling  theorem  are 
first  reviewed.  Using  the  idea  of  the  settling  time,  one  can  relate  the 
impulse  response  to  the  sampling  theorem.  Then  the  one  dimensional 
impulse  response  concept  is  extended  to  the  three  dimensional  space. 
Lewis  and  Bojarski's  [5,6]  work  in  the  area  is  also  briefly  contrasted. 

A  decision  rule  based  on  Petersen  and  Middleton's  [7]  definition  of 
sampling  efficiency  is  introduced  so  that  one  can  decide  on  the  more 
efficient  sampling  grid.  There  is  also  a  discussion  on  Petersen  and 
Middleton's  N  dimensional  sampling  theorem  and  its  different  forms  which 
depend  on  the  different  types  of  sampling  techniques.  Then  Mensa 
et  al.'s  [8]  polar  transformation  and  single  frequency  approach  to  two 
dimensional  Fourier  transform  is  rederived.  Their  approach  presents  a 
new  perspective  to  the  Nyquist  sampling  criteria. 

Some  practical  aspects  of  the  theory  in  Chapter  II  are  considered 
in  Chapter  III.  Even  though  the  discussion  is  in  one  dimension,  it  is 
also  applicable  to  two  and  three  dimensional  analyses.  The  real  signal 
requirement  on  Fourier  transform  is  rederived.  This  helps  to  clarify 
the  problem  of  non-zero  imaginary  parts  during  data  processing.  Using 
the  sampling  theorem  interpolation  scheme,  the  frequency  bandwidth's 


relationship  to  target  size  and  spatial  resolution  is  considered.  A 
common  problem  that  may  be  easily  overlooked  in  data  processing  is  the 
periodicity  of  the  Fourier  series  representation  to  solve  the  Fourier 
integral.  This  is  also  restated  in  the  last  part  of  Chapter  III. 

In  order  to  build  some  confidence  in  the  different  forms  of 
Petersen  and  Middleton's  sampling  criteria,  results  of  interpolating  one 
dimensional  impulse  responses  are  presented  in  Chapter  IV.  The 
interpolation  scheme  used  is  the  same  as  in  the  sampling  theorem,  except 
the  infinite  summation  is  a  finite  summation.  The  interpolated  results 
using  one  dimensional  sampled  data  appear  first  in  the  chapter  to 
preview  the  interpolation  via  two  dimensional  sampled  data.  Later,  an 
example  is  chosen  to  compare  the  efficiency  of  different  sampling  grids. 
Two  dimensional  impulse  responses  computed  using  discrete  two 
dimensional  Fourier  transform  are  compared  with  results  using  Mensa  et 
al.'s  approach  to  Fourier  transformation  in  Chapter  V.  The  potential  of 
using  the  spatial  impulse  response  for  target  imaging  is  explored  in  the 
last  part  of  the  report. 


CHAPTER  II 


THEORY 


Some  of  the  basic  concepts  on  impulse  response  and  one  dimensional 
sampling  theorem  are  individually  reviewed.  Then  the  two  concepts  are 
combined  and  extended  to  three  dimensional  space.  The  different  aspects 
of  N  dimensional  sampling  theorem  are  introduced.  Finally,  Mensa 
et  al.'s  approach  to  two  dimensional  Fourier  transform  and  sampling  is 
al so  di scussed. 

If  an  impulsive  electric  field  is  incident  on  a  target,  the 
normalized  far  zone  time  domain  scattering  will  be  the  impulse  response 
of  the  target  at  that  particular  aspect  angle  and  polarization.  The 
approach  to  be  used  herein  to  obtain  the  impulse  response  is  similar  to 
deconvolution.  First,  the  scattered  field  (a  complex  function  of 
frequency)  is  divided  by  the  spectrum  of  the  incident  wave.  This  result 
is  then  inverse  Fourier  transformed  to  produce  the  desired  impulse 


response. 


Let  f(t)  be  the  input  signal  in  time 

F(u>)  be  the  input  frequency  spectrum 
h ( t )  be  the  impulse  response  of  the  target 
H ( <*»)  be  the  frequency  response  of  the  target 
c(t)  be  the  output  signal  in  time 
C(w)  be  the  output  frequency  spectrum 


By  the  convolution  theorem  of  Fourier  transform  theory  [2] 


C(w)  =  F(«)H(u) 
C(w) 

<=>  H(u)  = 


Fourier  transformed. 


oo 
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h  (t )  -  2tt  / 

.OO 


C(O)) 


(2-1) 


(2-2) 


(2-3) 


According  to  Kennaugh  and  Moffatt  [3],  the  impulse  response  will 

r 

decay  exponentially  for  large  values  of  (t  -  ;r) 
where  t  -time 

r  -distance  between  observation  point  and  the  origin 
c  -speed  of  light 

Therefore,  one  can  define  a  settling  time  when  the  impulse  response  has 
its  amplitude  embedded  in  the  noise  level.  For  all  practical  purposes, 
the  settling  time  will  be  the  end  of  the  impulse  response.  Thus,  the 
impulse  response  is  a  time  limited  signal.  From  the  sampling  theorem 
[1,2]:  "A  time  limited  signal  can  be  reproduced  from  its  discrete 
frequency  values,  if  it  is  sampled  over  the  complete  frequency  domain 
using  the  Nyquist  rate."  As  a  consequence,  the  impulse  response  may  be 
reproduced  by  measuring  its  frequency  spectrum  at  the  Nyquist  sampling 
rate  (fs).  (i.e.,  f$  <  1/2T;  where  the  signal  is  limited  in  time  to 


f(t)  =  0 


It!  >  T 


then 


F(»)  = 


l 


nir  sin(coT  -  nn) 
f(T)  “3T  "-'"nir - 


n 


(2-4) 


where  the  samples  are  taken  at  «  =  nn/T,  but  n  is  taken  from  -00  to  ». 

Now,  let's  consider  the  concept  of  marching  in  time.  An  impulsive 
magnetic  field  incident  upon  a  solid  conducting  body,  sets  up  current 
J  on  the  surface.  As  a  result,  J  *  n  x  n  will  start  generating  a 

scattering  field  in  all  directions.  After  the  wave  passes  over  the 
target,  the  current  created  decays  exponentially.  The  scattered  field 
behaves  similarly.  At  one  aspect  angle,  the  time,  where  the  onset  of 
the  scattered  waveform  is  observed,  corresponding  to  the  time  required 
for  the  wave  to  reach  the  initial  scattering  point  on  the  target  and 
return  to  the  radar,  is  designated  as  the  initial  time  Ti#  The  time 
(Tf)  before  the  final  exponential  decay  occurrence  defines  the  end  of 
the  target. 

Let  x  be  the  length  of  the  object  along  the  line  of  sight 
at  one  aspect  angle 

At  be  the  time  difference  between  the  initial  time  and 
the  final  time  at  that  aspect  angle 
c  be  the  speed  of  light 


At  =  If  -  Ti 


then. 


2x  =  cAt 

x-c(Tf-Ti)  (2-b) 

“2 - 

Any  physical  object  has  finite  dimensions.  If  it  is  located  in  a 
rectangular  coordinate  system  (x:  x1,x2,x3),  then  it  can  be  said  to 

have  limited  dimensions  in  the  x  coordinates.  Namely,  the  object  is 
confined  by: 

va  b 

1  <  X1  <  X1 

a  b 

x2  <  x2  <  x2 

a  b 

x3  <  x3  <  x3 

.  a  a  a  b  b  b 
where  x^,  x2,  x3,  x^,  x2,  x3  are  some  real  constants. 

The  confinement  of  an  object  in  space  is  equivalent  in  saying  its 
spatial  impulse  response  is  time  limited,  as  the  impulse  response  has  a 
settling  time  for  every  aspect  angle.  Using  the  above  argument,  if  the 
impulse  response  measurement  is  obtained  at  all  aspect  angles,  then  an 
image  of  the  object  may  be  generated  from  the  data. 
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The  approach  taken  here  is  a  little  different  from  Lewi s-Bojarski ' s 
work  [5,6].  They  use  physical  optics  approximation  to  arrive  at  the 
formulation  for  imaging.  In  another  words,  they  do  not  use  the 
information  on  the  shadowed  side.  The  object  is  illuminated  in  every 
direction.  Information  on  the  4 it  solid  angle  is  required,  even  though 
measurement  on  the  complete  solid  frequency  sphere  is  difficult  to  be 
implemented.  With  the  help  of  the  N  dimensional  sampling  theorem,  which 
is  described  later,  the  implementation  becomes  more  practical.  If  one 
requires  information  over  a  finite  frequency  range,  then  the  infinite 
number  of  samples  over  the  4w  solid  angle  is  converted  to  a  finite 
number  of  samples.  The  N  dimensional  sampling  theorem  defines  a 
sufficient  sampling  criterion  to  characterize  space  limited  or  wave 
number  limited  signal.  Thus,  the  response  at  any  point  may  be 
interpolated  from  the  sampled  data  using  the  reconstruction  scheme 
defined  by  the  sampling  theorem. 

To  employ  the  sampling  theorem,  the  signal  must  be  limited  in  time 
or  frequency.  In  this  case,  the  object  is  limited  in  three  dimensions. 
The  transformed  space  will  be  the  wave  number  space.  Ideally,  the 
frequency  response  or  the  k-space  response  can  be  reproduced  by  sampling 
in  the  k-space  discretely  but  over  the  infinite  space.  To  obtain  the 
spatial  impulse  response,  a  three  dimensional  inverse  Fourier  transform 
is  performed  on  the  k-space  response. 


While  there  is  only  one  sampling  theorem,  there  are  different  forms 
of  the  sampling  that  satisfy  the  sampling  theorem.  In  one  dimension, 
the  different  forms  depend  on  how  the  time  limited  signal  or  the 
frequency  limited  signal  is  assumed  to  repeat  itself.  In  three 
dimensions,  the  forms  depend  on  how  the  target  is  placed  in  the 
reference  plane  (xj,  X2,  X3),  and  how  the  target’s  images  are  repeated 
in  the  three  dimensional  space.  Theoretically,  there  are  infinitely 
different  forms  of  the  sampling,  as  there  is  an  infinite  number  of 
different  target  shapes,  sizes  and  orientations.  One  would  prefer  a 
general  sampling  scheme  that  is  applicable  to  all,  or  at  least  most, 
objects  with  any  orientation.  This  is  where  the  canonical  containment 
cell  fits  in.  These  canonical  containment  cells  are  usually  of  simple 
geometries  so  that  a  wide  variety  of  targets  can  be  confined  by  their 
boundaries.  Examples  of  these  simple  geometries  are  sphere, 
parallelepiped,  ellipsoid,  finite  cylinder.  Thus  the  forms  of  the 
sampling  depend  on  the  choice  of  the  canonical  confinement  units  and 
how  the  unit's  images  are  repeated  in  the  three  dimensional  space. 

Two  simple  canonical  units  to  be  treated  in  this  report  are  the 
sphere  and  parallelepiped.  A  decision  rule  between  these  two  types  of 
units  will  be  discussed.  First,  the  concept  of  sampling  efficiency  will 
be  defined.  The  following  efficiency  formula  is  a  modified  version  of 
the  original  formula  defined  by  Petersen  and  Middleton  [7], 


CONSERVATIVE  ESTIMATE  OF  THE  VOLUME  OF  THE  OBJECT 
VOLUME  OF  THE  SMALLEST  SPHERE  ENCLOSING  THE  OBJECT 


=  Efficiency  for  isotropic  sampling 


CONSERVATIVE  ESTIMATE  OF  THE  VOLUME  OF  THE  OBJECT 
np  =  VOLUME  OF  THE  SMALLEST  PARALLELEPIPED  ENCLOSING  THE  OBJECT 

(2-7) 

=  Efficiency  for  parallelepipedic  sampling 


Rule: 


ns  >  np  use  spherical  confinement 


ns  <  np  use  parallelepiped  confinement 


(2-8) 


The  N  dimensional  sampling  theorem  obtained  by  Petersen  and 
Middleton  [7]  is:  "A  function  F(k)  whose  inverse  Fourier  transform  f(x) 
vanishes  over  all  but  a  finite  portion  in  x-space  can  be  everywhere 
reproduced  from  its  sampled  values  taken  over  a  lattice  of  points 


provi ded 


Ml  +  12“2  +  ...  +  ln“nh  U,  12,  ....  In  =  0,  ±  1,  ±  2,  ... 
that  the  vectors  {Uj}  are  small  enough  to  ensure  non-overlapping  of  the 
x-space  signal  f(x)  with  its  images  on  a  periodic  lattice  defined  by  the 
vectors  1^},  which  v^.Uj  =  2w6jj." 

Let's  consider  the  non-overlapping  condition  in  the  N  dimensional 
sampling  theorem.  The  requirement  is  'non-overlapping'  of  the  object 
cell  f(x)  with  its  images  on  a  periodic  lattice.  Thus,  the  periodic 
lattice  is  not  uniquely  defined.  Any  one  of  Figure  2-2,  2-3,  or  2-4  has 
a  valid  two  dimensional  periodicity.  Their  respective  periodicity  is 
defined  by  their  respective  N].t  ^2 } .  This  non-uniqueness  provides 
flexibility  on  the  choice  of  the  confinement  cell. 

However,  an  efficient  sampling  lattice  may  be  defined.  Pete,  sen 
and  Middleton  [7]:  "An  efficient  sampling  lattice  is  one  which  uses  a 
minimum  number  of  sampling  points  to  achieve  an  exact  reproduction  of  a 
space  limited  function."  In  another  words,  the  closest  packing  of  the 
object  cell  and  its  images  without  overlapping  will  be  efficient. 

^1,  v2»  v3^  will  be  changed  if  the  images  are  rotated  around  the 
object  cell;  hence,  the  sampling  lattice  is  still  not  fully  defined 
(Figures  2-2  and  2-4).  Nevertheless,  any  set  of  (vj,  V£>  V3}  defined  by 
the  above  criterion  will  have  the  same  efficiency,  or  the  same  number  of 
sampling  points. 


13 


The  application  of  the  sampling  theorem  for  the  reconstruction  of 
the  original  signal  needs  only  an  interpolation  formula,  provided  the 
non-overlapping  condition  is  met.  In  this  application,  if  one  has  the 
following: 

a)  A  signal  limited  in  x-space  and  its  images  are  specified 
by  vlf  v£  with  non-overlapping  condition  met. 

b)  The  sampling  lattice  in  k-space  is  defined  by  the  vector: 

Vi  +  12u2‘  where  *2  =  *  1.  *  2,  ±  3,  ... 

with  vi>u.  =  :  5^  =  Kronecker  delta 

or 

U  =  2nV-T 

where 

u  -  ca,|u'2] 

v  *  [5ji;2] 

-T  is  the  notation  for  the  transpose  of  the  matrix 
inverse  of  V 

c)  The  interpolation  formula 

then  one  can  reproduce  the  two  dimensional  impulse  response.  The 
procedures  are: 

1)  Sample  the  k-space  response  at  the  sampling  lattice  for 

FOlUl  +  1 2U2  ^  * 

2)  Interpolate  other  required  points,  if  necessary. 


The  interpolation  formula  taken  from  Petersen  and  Middleton  [7]. 


Ft*!,  Ic2)  =  I  l  F(l1iJ1  +  l2i72)G(k1k1  +  k2k2  -  ljUj  -  l2a2) 
1  ^=-00  1  2  =  -“ 

(2-9) 


where  k.  =  k^k^ 


3)  Inverse  Fourier  transform  the  k-space  response  to  obtain  the 
two  dimensional  impulse  response. 


f(xl*  x2}  =  4^  //  F(kpk2)eJ 


i(kl*l  + 


k  ^  x  ^  i 


dk^dk2 


(2-10) 


With  the  ease  of  calculation  in  mind,  the  periodicity  of  the  object 
cell  and  its  images  are  defined  as  in  Figure  2-5  for  this  report.  Its 
corresponding  sampling  lattice  is  shown  in  Figure  2-6. 

The  reconstruction  function  for  the  parallelogrammatic 
sampling  [7]: 


G((^,  u^) 


sin  iru),  sin  tto)0 

.  ( _ I)  ( _ i) 


TT 


IT  0^ 


(2-11) 


where 


is  in  the  direction  of  u^ 


id?  is  in  the  direction  of  U2 
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For  the  sphere  or  isotropic  confinement,  the  configuration  is 
adopted  from  the  concept  of  closest  packing  of  spheres  by  Coxeter  [9], 
The  configuration  chosen  for  this  report  is  as  shown  in  Figure  2-7.  The 
corresponding  sampling  lattice  is  shown  in  Figure  2-8. 

For  the  two  dimensional  case  [7,9]: 


*1 


=  R 


7 

-1 

7 


(2ir/R) 

”  2“ 
/3 

1 

✓5 

=  (2n/R) 

0 

L 

1 

(2-12) 


where  R  is  the  radius  of  the  isotropic  cell.  The  reconstruction 
function  [7]: 


1 

G(“l»  <*>2)  =  ~3  2  I  2.  x  t  2u>icos(R«i//7)cos(Rw2) 

R  ) 

-2usjC0S  (2Ra)j/ /3) 
-2y7w2Sin(Ru>i//T)sin(Rio2)  } 

(2-13)* 

where 

<*>1  is  in  the  direction  of  u^ 
is  in  the  direction  of  u^ 


(*  There  is  a  factor  missing  in  the  third  term  of  this  expression 
in  reference  [7],  See  Appendix  A  for  details  of  the  derivation). 


Since  this  report  only  deals  up  to  two  dimensional  sampling,  the  reader 
is  referred  to  Petersen  and  Middleton's  paper  on  the  three  dimensional 


reconstruction  function  G(io^,  u^,  io^). 

The  N  dimensional  sampling  theorem  gives  the  minimum  sampling 
lattice  or  criterion  which  will  sufficiently  define  the  k-space  signal. 
All  other  non  sampled  k-space  values  can  be  interpolated  via  an 
extension  of  Equation  2-9.  The  k-space  signal  can  be  inverse  Fourier 
transformed  into  the  spatial  domain  to  give  a  representation  of  the 
spatial  impulse  response.  In  another  words,  the  spatial  signal  is  also 
characterized  by  those  k-space  lattice  samples. 

If  one  is  interested  only  in  the  two  dimensional  impulse  response, 
then  Mensa  et  al .  [8]  presents  a  different  view  on  the  sampling 
criterion.  Unfortunately,  it  is  only  applicable  to  the  two  dimensions. 
First,  rectangular  to  polar  coordinate  transformation  is  applied  to  the 
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4n2 


<»  2  it 

/  /  rF(r,  ♦)ej2,'r,lcos^-9,dr  d« 


r=o  <j>=o 


(2-16) 


Then  the  double  integral  is  reduced  to  a  single  integral  by  considering 
only  a  particular  frequency  ring  (i.e.  S(r-<*>j)). 


“t 


2tt 


f i  (p»  9)  =  4tt2  /  F(o»1f  4> )e 
4>=o 


j2iro)^  pcos  ( 4> 


9) 

d<f> 


(2-17) 


Thus,  the  two  dimensional  Fourier  transformation  is  reduced  to  a 
convolution  type  integral.  F(kjt  k2)  is  the  frequency  response  of  the 
target  in  the  two  dimensional  k-space.  The  notion  of  S(r-uo^)  represents 
information  taken  only  with  one  frequency.  Subsequently,  the  two 
dimensional  impulse  response  is  obtainable  via  one  integration.  If 
information  from  other  frequency  rings  are  available,  then  superposition 
of  every  frequency  ring  response  in  the  spatial  domain  will  give  a  wide 
band  two  dimensional  impulse  response  representation. 

fT(xlf  x2)  *  l  f^P,  9)  (2-18) 

i 

where 

X1  =  PCOS0 
x2  =  psin© 

Naturally,  the  Nyquist  spacings  between  the  frequency  rings  and  between 
angular  samples  must  be  used,  before  the  total  two  dimensional  time 
response  obtained  can  be  considered  a  sufficient  representation  of  the 
true  two  dimensional  time  response. 
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The  angular  increment  which  satisfies  the  Nyquist  criterion  is 


given  by  Mensa  et  al .  [8]: 


X_ 
t  <  2D 


-1  X 

A0  <  <  sin  (2fr) 


TT 

v  <  4 


X  «  2D 

X  <  2D 

X  >  2D 


(2-19) 


where  D  =  maximum  dimension  of  the  object 

i 

X  =  wavelength  of  the  frequency  to  be  used 


(Assumption:  The  origin  of  the  x-space  coincides  with  the  middle  of  the 
object. ) 

The  author  would  like  to  propose  the  following  for  the  frequency 
i ncrement : 


c 


Af  <  2(1  '+  k)b 


(2-20) 


where  D  =  maximum  dimension  of  the  object 
c  =  speed  of  light 
K  =  some  safety  factor 


The  aspect  angle  which  has  the  longest  dimension  of  the  object  is 
assumed  to  have  the  longest  settling  time  in  its  impulse  response.  All 
other  aspect  angles  require  Nyquist  frequency  increments  larger  than  or 
equal  to  this  aspect  angle.  The  value  of  the  safety  factor  is  the  best 
estimation  achievable  by  other  means.  The  reason  for  the  safety  factor 
is  not  all  impulse  response  signals  are  limited  to  the  length  of  the 
object  at  any  aspect  angle.  Furthermore,  in  the  GTU  sense,  some  of  the 


multiple  scattering  effects  may  be  eliminated  or  included  by  employing 
smaller  or  larger  value  of  the  safety  factor.  This  is  equivalent  to 
truncation  of  the  signal  in  time  during  measurement. 


CHAPTER  III 


PRACTICAL  CONSIDERATIONS 

In  this  chapter,  the  practical  aspects  of  the  previously  described 
theory  are  presented.  Since  negative  frequencies  cannot  be  physically 
generated,  an  assumption  on  the  negative  frequency  response  must  be 
made.  First,  the  assumption  of  a  real  measured  signal  is  discussed. 

The  frequency  limits  are  then  considered  in  relation  to  the  object's 
size,  its  impulse  response  and  the  data  processing  requirement.  For 
clarity,  the  practical  aspects  are  discussed  in  one  dimension. 

Extension  to  the  higher  dimensions  can  be  easily  accomplished. 

In  general,  a  true  impluse  is  hard  to  generate;  instead,  a  Gaussian 
pulse  is  often  used.  There  are  also  times  when  good  narrow  Gaussian 
pulses  are  not  readily  available.  In  these  cases,  the  approach  of 
frequency  sweeping  may  be  used.  The  bandwidths  of  most  oscillators, 
waveguide  components,  transmitting  and  receiving  antennas  are  limited. 

In  essence,  the  frequency  band  can  only  be  swept  from  ^  to  “h  (Figure 
3-1).  To  overcome  part  of  the  problem,  let's  consider  the  one 
dimensional  sampling  theorem  [1,  2]: 


If  F(<o)  is  even,  then  the  integrand  in  the  second  integral  is  an 
odd  function  of  w.  it  follows  that  the  second  integral  is  zero. 
Therefore, 


1 

f(t)  *  2n  /  F(to)cos(u>t)  du>  (3-4) 

-00 

=  real  function,  if  F(o>)  is  even 

(Note:  If  F(w)  is  complex,  then  Re[F(w)]  =  Re[F(-to)]  and 
Im[F(<o)]  =  -Im[F(-u)]  are  the  conditions  for  f(t) 
to  be  real ). 

In  this  discussion,  the  object  is  real.  It  follows  that  f(t)  is 
also  real.  Now,  information  from  and  \  **  10  ^  :*V|  i  s 

available.  (Figure  3-2) 

Next,  the  Rayleigh  Law  is  used  to  determine  the  scattered  field  at 
d.c.  or  zero  frequency.  From  Kennaugh  and  Cosgriff  [3]:  "As  the  source 
frequency  tends  to  zero,  all  finite  scatterers  follow  the  Rayleigh  Law, 
giving  a  scattered  field  intensity  which  diminishes  as  the  square  of 
frequency."  The  scattered  field  intensity  at  d.c.  will  be  zero. 

If 


IT 

T  > 


then 


Fn  is  known  for  -N  <  n  <  N 


-HI, 


Figure  3-2.  Frequency  information  available  after  assumption  is  made 
about  the  negative  frequency  samples 
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where 


N  =  I  [—  3  (3-5) 

I[x]=truncation  of  x's  value  after  the  decimal  point 

IT 

That  is,  if  the  first  sampled  location  (y)  is  higher  than  or  equal  to 
UL,  then  information  is  available  from  to  because  one  can 
interpolate  the  in  between  data  points  using  the  sampling  theorem. 

Let 

*s  be  the  object  size 
K  be  some  safety  factor 

then  the  settling  time  Ts  is 

T  =  2(1+K)xs  with  T  =  T  /2 

IT- 

One  must  have  the  condition: 

IT 

T  >  \ 

2irc 

2( 1+K)xs  >  \ 

1TC  2lTC 

<=>  (1+K)xs  >  XL 

<=>  \  >  2(1+K)x$  (3-6) 

This  puts  a  limit  on  the  lowest  frequency  usable,  or  the  largest  object 
size  assumed  by  a  given  W|_. 


33 


The  span  effectively  determines  the  resolution  of  the 

impulse  response  signature.  Let  xr  be  the  desired  resolution  on  the 
object,  then  the  impulse  response  resolution  is 

2xr 
tr  =  c 

2  ire 

<=>  w  =  9— 

r  r 

One  has  the  condition: 

2itc 

*H  >  2^:  =  “r 

2irc  ire 

<=>  "Th  >  xF 

<=>  AH  <  2*r  (3-7) 

Conditions  3-6  and  3-7  will  help  to  decide  how  wide  a  frequency  band  may 
be  used.  If  the  bandwidth  is  wide  enough,  then  the  impulse  response  can 
be  generated  to  a  very  good  approximation. 

A  common  problem  during  implementation  may  involve  the  inverse 
Fourier  transform  on  the  reconstructed  frequency  spectrum.  This 
waveform  may  or  may  not  be  a  well  defined  function  which  can  be  inverse 
Fourier  transformed  into  a  closed  form  solution.  The  common  approach 
would  be  to  approximate  the  integration  using  a  summation  on  a  digital 
computer.  In  effect,  this  approach  will  be  a  Fourier  series 
representation  which  requires  the  time  or  frequency  waveform  to  be 
periodic.  Consequently,  there  are  2  more  limitations: 


1)  The  period  (Tc)  used  in  the  digital  computation  must  be  greater 
than  two  times  the  settling  time  (Ts)  of  the  impulse  response. 

( i • e . ,  Tc  >  2Tj) 

2)  The  period  (o^-)  used  in  the  digital  computation  must  be  greater 
than  two  times  the  highest  frequency  (oj^)  swept. 

(i  .e.,  >  2ufl) 

Fortunately,  a  standard  IBM  subroutine  FFT  (Fast  Fourier  Transform) 
package  is  available  to  do  the  required  Fourier  analyses. 

Furthermore,  one  dimensional  impulse  response  requires  a  two 
dimensional  plot.  The  two  axis  quantities  are  amplitude  and  time.  Two 
dimensional  impulse  response  requires  a  three  dimensional  plot.  The 
three  axis  quantities  are  the  amplitude  and  the  plane  axes.  Should 
anyone  consider  three  dimensional  impulse  response,  one  requires  a  plot 
in  four  dimensional  space.  Consequently,  this  report  will  only  deal 
with  impulse  responses  up  to  two  dimensions.  With  the  knowledge  in 
one's  mind  that  the  spatial  impulse  response  can  easily  be  obtained  by 
an  extension  of  this  two  dimensional  approach  when  there  is  an 
appropriate  representation. 
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CHAPTER  IV 


ONE  DIMENSIONAL  IMPULSE  RESPONSES 


In  this  chapter,  the  interpolation  of  one  dimensional  impulse 

i 

1  responses  is  presented;  first,  the  result  of  interpolation  using  one 

dimensional  data;  then,  using  two  dimensional  data.  The  object  is  a  six 
inch  diameter  metallic  sphere.  This  object  choice  is  because  the  Mie 

! 

1  solution  in  frequency  domain  is  readily  available.  In  this  first 

section,  the  Mie  solution  frequency  data  are  sampled  at  different  rates 
and  interpolated  either  using  a  straight  line  or  a  sine  reconstruction 
function. 


F(»)  = 


N 

l 

n=-N 


,nir 

F(T) 


si  n  ( o>T  -  n^) 
wT  -  n* 


where 


(4-1) 


,jJHT 

N  =  1  (  IT  ) 

I ( x )  =  truncation  of  x's  value  after  the  decimal  point 
UIH  =  highest  frequency  used 


[Note:  This  is  Equation  (3-1)  except  the  summation  is  from 
-N  to  +  N], 


A  cosine  tapering  weighting  function  [10]:  (Figure  4-1)  (for  the 
convenience  of  the  reader,  all  tables  and  figures  of  Chapter  IV  are 
grouped  together  at  the  end  of  the  chapter) 


W(n), 

L 


irn 

=  0.5[1+cos(nj)] 
=  0 


I  n  |  <  Ni 

1 n  t  >  Ni 


(4-2) 


is  multiplied  to  the  interpolated  frequency  data.  The  purpose  of  the 
weighting  or  filtering  is  to  reduce  the  effect  of  the  Gibb's 
phenomenon  [2],  The  resulting  data  are  inverse  Fourier  transformed  into 
the  time  domain  discretely  to  give  the  impulse  response  in  time. 

Figure  4-2  represents  a  time  domain  impulse  response  plot  obtained 
from  the  Mie  solution  for  a  six  inch  diameter  metallic  sphere  using  the 
frequency  spectrum  from  0  to  12  Ghz  with  a  cosine  tapering  filter.  The 
imperfect  specular  impulse  and  the  ripples  around  it  at  the  start  of  the 
response  are  caused  by  1)  finite  bandwidth,  and  2)  Gibbs'  phenomenon. 

Figure  4-2  is  considered  to  be  the  standard  for  comparison  with  the 
other  responses.  Its  frequency  samples  are  taken  every  60  Mhz  and 
linked  together  by  straight  lines.  Figure  4-3  has  frequency  samples 
taken  every  0.125  Ghz  over  the  spectrum  of  0  to  12  Ghz  and  interpolated 
the  in  between  points  using  Equation  (4-1).  The  differences  of  the 
impulse  responses  in  this  chapter  from  Figure  4-2  (the  'exact'  solution) 
are  shown  in  figures  designated  with  their  respective  figure  number  plus 
an  'a'  attached.  For  example,  to  obtain  Figure  4-2  from  Figure  4-3,  one 
adds  Figure  4-3a  to  Figure  4-3. 
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Figures  4-4,  4-5,  and  4-6  are  similar  to  Figure  4-3,  except 
frequency  samples  are  taken  every  0.25  Ghz,  0.5  Ghz,  and  0.76  Ghz 
respectively.  Figure  4-5  is  still  recognizable  to  be  Figure  4-2  but 
Figure  4-6  is  not.  Figure  4-6  does  not  have  similar  behavior  because  it 
employs  a  frequency  sampling  rate  below  the  Nyquist  rate.  From  Figure 
4-2,  one  can  estimate  the  settling  time  of  the  impulse  response  of  the 
sphere  to  be  about  1.5  nsec.  Equally  well,  one  could  use  a  longer 
settling  time  depending  on  one's  assumption  of  the  noise  amplitude, 
i  .e. , 

T  =  0.75E-9S 

1  1 
fs  >  2T  -  1.5E-9S 

=  0.66  Ghz 

<  0.75  Ghz  (Figure  4-6) 

Figures  4-2  to  4-5  have  sampling  rate  better  than  the  Nyquist  rate. 
If  the  sampling  rate  is  high  enough,  then  one  may  use  straight  line 
interpolation  (Figure  4-2)  instead  of  the  sine  reconstruction  function 
to  save  computer  time  on  interpolation.  On  the  other  hand,  if  samples 
are  scarce  but  still  satisfy  the  Nyquist  criterion,  then  the  sine 
reconstruction  function  (Equation  (4-1))  is  preferred  for  more  pleasing 
results.  (Figures  4-3,  4-4,  4-5) 

Now,  the  reconstruction  of  the  one  dimensional  impulse  responses 
from  data  obtained  on  two  dimensional  isotropic  and  cubic  sampling 
lattice  is  presented. 

A  A 

F(k^,  kg)  =  I  I  F(l^uj  +  ^2u2^^1*c1  +  ^2^2  ”  ^lul  ”  ^2U2^ 

]1  ]2  (4-3) 
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where 


•  -  12  are  summed  over  1^  and  12  which  satisfy 

kl  <  1 1 lu 1  +  1 2U2 1  <  kH 

kL  and  kH  define  the  frequency  bandwidth  used. 

(Note:  This  is  Equation  (2-9)  with  finite  summation.) 


Samples  are  taken  out  of  the  Mie's  frequency  solution  on  the  isotropic 
sampling  lattice  defined  by  u^  and  U2  in  Equation  (2-12),  and  the  cubic 
sampling  lattice  defined  by  : 

U1  =  2*  [J]  u2  =  2 it  [°]  (4-4) 

i  .e. , 


samples  are  taken  over  1^  +  l2u2  for  lj,  12  =  0,  ±  1,  ±  2,  ... 

(4-5) 

This  sampling  lattice  is  generated  by  the  program  PTGRID  (see  Appendix 
B).  The  frequency  response  for  an  aspect  angle  is  interpolated  using 
Equation  (4-3)  in  conjunction  with  either  Equation  (2-13)  for  isotropic 
sampling  or  Equation  (2-11)  for  cubic  sampling.  This  interpolation  work 
is  done  by  the  program  INTERPOL  (see  Appendix  B).  The  resulting 
frequency  response  is  again  cosine  tapering  low  pass  filtered  to  reduce 
the  Gibb's  phenomenon.  After  the  filtering,  the  frequency  spectrum  is 
inverse  Fourier  transformed  discretely  into  the  time  domain  to  give  an 
impulse  response  picture  for  the  six  inch  metallic  sphere.  The 


filtering  and  the  discrete  inverse  Fourier  transform  are  functions  of 
FTRAN  [11]  (see  Appendix  B). 

Sampling  at  every  0.6  Ghz  for  a  six  inch  diameter  metallic  sphere, 
is  equivalent  to  assuming  a  signal  having  four  times  the  diameter  of  the 
sphere  in  one  dimensional  space.  Therefore,  the  two  dimensional 
confinement  cell  is  assumed  to  contain  the  sphere  and  has  a  guard  band 
of  1.6  times  the  maximum  dimension  of  the  sphere  surrounding  the 
sphere.  The  safety  factor  thus  chosen  is  1  (i.e.,  2(1+K)=4  <=>  K=l). 

The  frequency  range  sampled  is  0  to  12  Ghz. 

Figures  4-7,  4-8,  4-9,  and  4-10  are  one  dimensional  impulse 
responses  reconstructed  from  two  dimensional  isotropic  sampling  data, 
for  aspect  angles  of  0,  0.719,  1.438,  30  degree  respectively.  The  CPO 
time  taken  for  interpolating  each  waveform  is  about  6  minutes  for  200 
plotting  points.  Figures  4-11,  4-12,  4-13,  4-14  are  reconstructed  from 
cubic  sampling  data,  for  aspect  angles  of  0,  1.193,  2.386,  45  degrees. 
The  CPO  time  taken  for  these  waveforms  is  about  2.5  minutes  for  200 
plotting  points.  Less  time  in  interpolation  for  cubic  sampled  data  is 
probably  due  to  the  simplicity  of  the  reconstruction  function  (Equation 
(2-11)  versus  (2-13)).  These  angular  choices  are  arbitrarily  chosen. 

One  should  note  the  close  resemblance  of  all  these  figures  (4-7  to  4-14) 
with  Figure  4-5.  Next,  let's  consider  the  case  of  more  samples  taken. 
Effectively,  the  safety  factor  is  changed  from  one  to  three  but  the 
frequency  range  remains  the  same.  Figures  4-15,  4-16,  and  4-17  are 
reconstructed  one  dimensional  impulse  responses  using  isotropic  samples 
for  aspect  angles  of  0,  0.352,  30  degrees  respectively.  Again 
discrepancy  is  not  high  (Figures  4-15a,  4-16a,  and  4- 1 7 a ) .  Since  the 
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density  of  samples  taken  is  finer,  or  more  samples  participated  in  the 
interpolation,  the  CPU  time  has  incresed  to  about  20  minutes  for  200 
plotting  points. 

If  both  the  isotropic  and  cubic  sampling  can  perform  competitively, 
how  does  one  decide  on  which  sampling  grid?  Tne  answer  lies  in  the 
efficiency  definition  defined  previously.  However,  the  area  is 
considered  in  a  plane  instead  of  the  volume  in  a  three  dimensional 
space. 

i.e.,  one  modifies  Equations  (2-b)  and  (2-7)  to 

AREA  OF  A  CROSS-SECTION  ON  THE  OBJECT 

Efficiency  -  AKEA  UF  THE  SAME  CROSS-SECTION  ON  THE  TWO  (4-6) 

DIMENSIONAL  ENCLOSURE 

but  the  decision  rule:  Equation  (2-8),  remains  the  same. 

Efficiency,  as  mentioned  before,  is  defined  as  a  minimum  sampling 
requirement.  Three  objects:  a  six  inch  diameter  sphere,  a  3/7  inch 
cube  and  a  sphere  cap  cylinder  are  shown  in  Figure  4-18  on  their  major 
axis  cross-section.  Their  respective  cross-sectional  areas;  closest 
circular,  squared,  rectangular  enclosure  cross-sectional  areas;  and 
efficiencies  are  tabulated  in  Table  4-1. 

The  definition  of  Equation  (4-6)  is  used  to  locate  the  sampling 
lattice.  The  number  of  these  locations  over  a  frequency  range  is  summed 
to  give  the  minimum  number  of  sampling,  defined  by  the  sampling  theorem, 
to  sufficiently  characterize  the  spatial  impluse  response  in  that 
frequency  range.  This  work  is  done  by  the  program  PTGR1D  (see  Appendix 
B).  The  numbers  are  tabulated  in  Table  4-2  and  plotted  in  Figures  4-19, 


41 


4-20,  and  4-21.  A  safety  factor  of  one  is  used  in  the  computation.  The 
sizes  of  the  objects  are  chosen  such  that  the  circular  enclosure  is  the 
same  for  all  three  objects  for  easy  comparison  in  the  graphs. 

As  the  frequency  range  becomes  larger  and  larger,  the  significance 
of  the  efficient  sampling  grid  becomes  more  and  more  important.  Let's 
take  the  example  of  the  sphere  cap  cylinder  (Figure  4-19).  The  use  of 
the  rectangular  enclosure  will  provide  an  efficiency  of  0.96.  The 
number  of  samples  required  over  0  to  12  Ghz  is  696.  The  use  of  the 
squared  enclosure  can  only  give  an  efficiency  of  0.38.  Tne  number  of 
samples  required  is  2.6  times  that  of  the  rectangular  enclosure.  The 
use  of  the  circular  enclosure  has  an  efficiency  of  0.45.  The  number  of 
samples  is  about  2.3  times  that  of  the  rectangular  enclosure.  As  the 
frequency  range  is  expanded,  saving  in  measurement  time  by  the  proper 
choice  of  the  enclosure  becomes  very  substantial. 


Object 

SPHERE  ] 

CUBE 

SPHERE  CAP  CYLINDER 

Area  on 
major  axis 

182.4 

l 

116.1 

82.2 

Area  of  closest 
ci rcular  enclosure 

182.4 

i 

182.4 

182.4 

Area  of  closest 
squared  enclosure 

232.3 

! 

lib.  1 

214.7 

Area  of  closest 
rectangular  enclosure 

232.3 

116.1 

86.9 

Efficiency  of 
ci rcular  enclosure 

1 

0.64 

0.46 

Efficiency  of 
squared  enclosure 

0.79 

1 

0.38 

Efficiency  of 
rectangular  enclosure 


0.79 


1 


U.9fa 


TABLE  4-2 

THE  NUMBER  UF  SAMPLING  POINTS  FOR  3  TYPES  OF  SAMPLING 
OU  3  OBJECTS  OVER  VARIOUS  FREQUENCY  RANGES 


Object 

SPHERE 

CUBE 

SPHERE  CAP  CYLINDER 

Type  of  enclosure 

ci rcular 

squared 

squared 

rectangul ar 

squared 

Frequency  ranges 
from  0  up  to 
(Ghz) 

1 

12 

12 

8 

2 

8 

2 

42 

48 

24 

20 

44 

3 

9b 

120 

60 

46 

108 

4 

186 

212 

100 

8U 

192 

S 

282 

324 

160 

lib 

292 

b 

408 

472 

240 

174 

436 

7 

bb8 

640 

324 

592 

8 

720 

828 

420 

768 

9 

912 

1048 

516 

384 

972 

10 

1134 

13U4 

666 

476 

1200 

11 

1368 

1564 

776 

582 

1456 

12 

1626 

1876 

940 

696 

1740 

4-19 

- 

- 

4-19 

4-19 

Plotted  in 

4-20 

4-20 

- 

- 

Fi gure 

4-21 

“ 

4-21 

“ 

Notation : 

X 

0 

0 

a 

0  .1 

I) 


NTH  HARMONICS 


Cosine  tapering  weighting 
(Equation  (4-2 )  with  =  16) 


4-3.  Impulse  response  of  a  6"  metallic  sphere  with  frequency 
samples  taken  every  125  Mhz  over  the  range  of  0  to  12  Ghz 
and  interpolated  using  Equation  (4-1) 
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Figure  4-4.  Similar  description  as  Figure  4-3,  except  frequency  samples 
are  taken  every  2S0  Mhz 
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igure  4-b.  Similar  description  as  Figure  4-3,  except  frequency  samples 
are  taken  every  500  Mhz 
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.  Impulse  response  of  a  6"  'tietallic  sphere*  at  0°  aspect  angle 
with  the  1-D  frequency  response  interpolated  from  2-1) 
isotropic  lattice  samples  taken  with  a  safety  factor  of  1 
over  the  range  of  0  to  12  Ghz  using  tquation  (4-3)  and 
and  (2-13) 
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Figure  4-7a,  Error  of  Figure  4-7  from  Figure  4-2 
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Figure  4-10.  Similar  description  as  Figure  4-7,  except  the  aspect  angle 
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re  4-11.  Impulse  response  of  a  6"  metallic  sphere  at  0°  aspect 

angle  with  the  1-0  frequency  response  interpolated  from 
2-0  cubic  lattice  samples  taken  with  a  safety  factor 
of  1  over  the  range  of  0  to  12  Ghz  using  Equation 
(4-3)  and  (2-11) 
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4-15.  Impulse  response  of  a  6"  metallic  sphere  at  0°  aspect 

angle  with  the  1-0  frequency  response  interpolated  from 
2-D  isotropic  lattice  samples  taken  with  a  safety  factor 
of  3  over  the  range  of  0  to  12  Ghz  using  Equation  (4-3) 
and  (2-13) 


TIME  IN  NBNOSECS 


riptlon  as  Figure  4-15,  except  the  aspect 
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Error  of  Figure  4-17  fron  Figure  4-2 
(Magnification  =  lOx) 
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Figure  4-iy.  The  number  of  samples  specified  by  the  different  types  of 
sampling  lattices  on  a  sphere  cap  cylinder  described  in 
Figure  4-18 
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CHAPTER  V 


TWO  DIMENSIONAL  IMPULSE  RESPONSES 

Thus  far,  the  discussion  has  focused  only  on  the  number  of  samples 
required.  The  smaller  the  number  of  samples,  the  less  measurement  time 
is  required.  However,  if  the  two  dimensional  impulse  response  is  of 
interest,  then  the  transform  method  must  also  be  considered.  This 
chapter  will  discuss  the  potential  time  consumption  problem  of 
multi -dimensional  Fourier  transform  plus  some  possible  solutions.  Image 
reconstruction  using  the  spatial  impulse  response  is  also  considered. 

While  an  isotropic  enclosure  may  be  more  efficient  to  enclose  a 
sphere  than  a  cube,  it  is  not  as  easy  to  do  three  dimensional  discrete 
Fourier  analyses  as  the  cubic  enclosure.  Most  of  the  discrete  Fourier 
transform  techniques  are  developed  to  fit  equally  spaced  data.  In 
another  words,  programs  are  written  to  perform  readily  on  the  cubic 
sampling  lattice.  Any  other  sampling  grid  data  must  be  interpolated  to 
the  cubic  grid  either  before  or  after  the  discrete  Fourier  transform 
step;  otherwise,  the  proper  representation  cannot  be  achieved. 

Sometimes,  interpolation  may  also  be  desired  for  the  cubic 

sampling  grid  data,  as  in  the  case  of  higher  resolution  requirement  than 


measured.  This  type  of  interpolation  reqjires  a  lot  of  computer  time  as 
pointed  out  in  Mensa  et  al.'s  paper  [8], 

Let's  consider  an  example  to  estimate  the  time  involved. 

From  Chapter  IV: 

1)  Interpolation  for  2U0  plotting  points  using  1626 
isotropic  samples  is  about  6  mins. 

2)  Interpolation  for  20U  plotting  points  using  1876 
cubic  samples  is  about  2.5  mins. 

3)  The  number  of  samples  required  on  a  sphere  cap 
cylinder  for  circular,  rectangular  and  squared 
enclosures  are  1626,  696,  174U  respectively. 

(Table  4-2) 

Compact  range  measurement  facility  at  0.  S.  U.  per  Walton  [12 J: 

The  measurement  system  response  time  is  about: 

1)  0.2s/data  point,  if  frequency  scan  is  used. 

2)  1  min/frequency  ring,  if  angular  scan  is  used. 

The  time  estimation  for  each  type  of  the  sampling  measurement  on  the 
sphere  cap  cylinder  and  two  dimensional  interpolation  is  presented  in 
Table  5-1.  (For  the  convenience  of  the  reader,  all  tables  and  figures 
of  Chapter  V  are  grouped  together  at  the  end  of  the  chapter.)  The 
interpolation  time  for  the  rectangular  enclosure  data  may  also  be 
considered  zero.  By  remembering  a  scale  factor,  the  data  can  be 
processed  as  in  a  squared  grid.  If  the  proper  signal  representation  or 
a  finer  resolution  is  required,  the  interpolation  step  is  still 


unavoidable.  Therefore,  the  numbers  in  the  table  do  give  a  fair 
comparison,  as  all  data  are  brought  to  the  same  level  -  squared  grid 
representation.  The  interpolation  step  plays  an  important  role  on  the 
decision  of  which  type  of  sampling  to  use.  The  saving  in  measurement 
time  is  sometimes  balanced  out  by  the  data  processing  time. 

The  Nyquist  criteria  (Equations  (2-19),  (2-20))  plus  the  polar 
transformation  (Equation  (2-15))  described  earlier  give  new  perspective 
for  the  efficient  but  non-cubic  enclosure.  Now  one  only  needs  to 
interpolate  for  a  sufficient  number  of  frequency  rings.  This  work  is 
done  by  program  INTERPOL  (see  Appendix  B).  Each  frequency  ring  is 
transformed  to  the  spatial  domain  individually  and  summed  together  to 
get  the  total  time  response.  The  former  is  the  work  of  program  INTEGFFT 
(see  Appendix  B);  the  latter  is  the  job  of  program  SUM3D  (see  Appendix 
B).  The  interpolation  performed  this  way  may  require  as  much  time  as 
the  interpolation  onto  the  squared  grid.  It  nevertheless  gives  an 
alternative  way  to  the  solution  of  the  problem.  Now,  another  question 
may  be  raised:  Why  perform  interpolation  if  it  is  so  time  consuming? 

The  interpolation  may  be  avoided,  if  the  two  dimensional  time 
response  is  the  only  interest.  All  one  has  to  do  is  to  measure  the 
frequency  rings  and  then  process  the  data  as  described  before.  However, 
if  one  desires  the  impulse  response  of  a  target  at  a  particular  aspect 
angle  or  the  response  over  a  frequency  ring  other  than  those  measured, 
one  is  required  to  develop  a  different  interpolation  scheme  than  the  one 
described  in  the  theory  section.  A  possible  way  to  reduce  the 


interpolation  time  is  to  derive  a  faster  interpolation  scheme  or  a 
general  multi -dimensional  Fourier  transform  technique  which  operates  on 
any  grid.  A  parallel  processor  which  uses  optics  may  be  another 
possibility  in  terms  of  hardware.  A  lens  system  has  a  Fourier 
transform  relationship  between  the  source  and  the  image  [13], 

The  Mie  solution  for  a  sphere  is  a  very  good  data  example;  except 
its  backscattered  response  is  isotropic.  A  proper  choice  of  test  object 
must  have  non-isotropic  property.  The  sphere  shifted  off  the  centre  of 
the  plane  is  one  possibility,  but  it  only  has  variation  in  the  phase 
term  and  not  in  the  magnitude  term.  Since  most  of  the  other  exact 
solutions  are  not  readily  available,  a  first  order  UTD  solution  for  a 
finite  circular  cylinder  [14]  is  used;  with  the  caution  that  UTD  gives  a 
valid  approximation  to  the  exact  solution  at  high  frequency. 

The  size  of  the  circular  cylinder  is  chosen  to  be  six  inches  in 
length  and  three  inches  in  diameter.  Having  chosen  the  size  of  the 
cylinder,  one  has  defined  the  low  frequency  limit  of  this  UTD  model  to 
about  2  Ghz.  Since  the  step  and  ramp  response  of  an  object  reduce  the 
need  for  the  high  frequency  information,  only  the  impulse  response  of 
this  cylinder  solution  is  considered  here. 

Figure  5-1  is  the  two  dimensional  impulse  response  (Equation  (2-9)) 
for  the  above  circular  cylinder.  The  frequency  samples  are  first 
generated  on  the  grid  points  defined  by  Equation  (4-5)  over  only  180°. 

In  this  case,  the  spatial  impulse  response  is  assumed  to  settle  after 
the  wave  has  passed  over  four  times  the  length  of  the  object  at  every 
aspect  angle.  The  rectangle  so  designated  has  dimensions:  24"  in 


length  and  12"  in  diameter.  The  periodic  lattice  is  chosen  as  in  Figure 
2-3.  Consequently, 

\  =  24"  [J]  v2  =  12"  [*] 

The  guard  band  is  9"  on  each  side  of  the  cap  of  the  cylinder  and  4.5" 
on  the  circular  surface.  Then, 


2tt 

U1  =  24" 


ll_  rO, 

u2  =  12"  L 1  * 


Since, 


1  lui  =  hlullkl 

/V 

l2u2  =  1 2  ^  u2 i k2 

where 

1 x,  12  =  0,  ±  1,  ±  2,  ±  3,  ... 

Therefore,  the  sampling  lattice  is  defined  by: 

A  A 

1  lu  1  +  12U1  =  1 1 1 U1 1  kl  +  ^  2 1  u2  i  k2 
<=>  points  on  the  k-plane: 


kl  =  hlull  and  k2  =  1 2 1 u2 1 


8b 


where 


6 

i 

) 


» 


» 


1 1,  12  =  U,  ±  1,  ±  2,  ±  3,  ... 

Then  the  sampled  data  are  Fourier  transformed  discretely  into  the 
spatial  domain.  The  discrete  two  dimensional  Fourier  transform  is 
performed  by  an  IBM  subroutine  HARM  (see  Appendix  8).  The  frequency 
range  used  is  2  to  b  Ghz.  The  safety  factor  of  one  is  used.  The 
vertical  polarization  data  are  taken  from  the  UTD  solution  on  a  major 
axis  cut.  Figure  b-la  is  the  contour  plot  of  Figure  b-1.  Figure  b-2  is 
also  a  two  dimensional  impulse  response  for  a  cylinder  solution  that  has 
the  same  previous  description.  The  differences  are  the  technique  of  the 
Fourier  transformation  and  the  sampling  locations  tn  the  frequency 
plane.  It  employs  Mensa  et  al.'s  method  on  6  different  frequency  rings: 
2.5,  3,  3.5,  4,  4.5,  5  Ghz.  (Equation  (2-18))  The  samples  are  taken  so 
that  Equations  (2-19)  and  (2-20)  are  satisfied.  Again  the  samples  are 
taken  over  180°.  Figure  5-2a  is  the  contour  plot  of  Figure  b-2. 

Comparing  Figures  5-1  and  b-2,  one  can  see  many  similarities.  The 
differences  can  be  deduced  by  recalling  their  respective  generating 
methods.  Figure  5-1  has  information  scattered  all  over  the  frequency 
range;  while.  Figure  b-2  has  vaules  only  over  those  frequency  rings 
mentioned  before.  There  is  also  the  processing  error  involved. 

One  interesting  thing  is  to  be  noted  in  Figure  b-1,  or  b-2.  If  one 
records  just  the  highest  points  within  its  neighborhood,  one  can  trace 
out  a  rectangle.  This  may  be  easier  to  see  on  a  contour  plot  (Figure 
b-la,  5-2a).  This  looks  like  one  cross-section  of  the  target.  If  2" 
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solid  angle  information  is  available,  then  an  image  of  the  target  can 
indeed  be  produced.  Of  course,  the  resolution  is  still  governed  by  the 
highest  frequency  used.  In  this  case  the  bandwidth  used  is  quite 
narrow;  hence,  the  resulting  resolution  is  not  high. 

Figure  5-3  is  generated  similarly  as  Figure  5-2.  It  is  generated 
using  Mensa  et  al.'s  method  on  16  different  frequency  rings:  2.5,  3, 
3.5,  4,  4.5,  5,  5.5,  6,  6.5,  7,  7.5,  8,  8.5,  9,  9.5,  10  Ghz.  Figure 
5-3a  is  the  contour  plot  of  Figure  5-3.  The  dimensions  of  the  rectangle 
defined  by  the  high  peaks  do  correspond  to  the  major  axis  cross-section 
of  the  circular  cylinder  described  earlier.  One  may  wonder  how  a 
rectangle  is  concluded  from  Figure  5-3a.  There  are  a  few  clues 
available.  The  Tj's  are  quite  obvious  from  the  high  peaks  on  the 
illuminated  side.  The  final  time  (Tf)  of  the  one  dimensional  impulse 
response,  as  one  recalls,  is  defined  by  the  beginning  of  the  exponential 
decay  of  the  signal.  The  final  peaks  at  the  coordinates  (5,3)  and 
(5,-3)  are  indeed  higher  than  any  point  x^  >  5.  The  ridges  and  valleys 
on  the  shadow  side  (xi  >  5)  of  the  cylinder  is  fairly  straight. 

Let's  consider  the  case  where  this  imaging  theory  converges  to 
Lewi s-Bojarski ' s  work.  Figure  5-4  is  the  same  as  Figure  5-3  except 
data  are  taken  over  360°,  or  illuminated  in  every  direction  on  the  two 
dimensional  plane.  Figure  5-4a  is  the  contour  plot  of  Figure  5-4.  As 
expected,  the  peaks  (or  T^'s)  in  the  figure  trace  out  a  rectangle  which 
is  the  major  axis  cut  of  the  circular  cylinder.  One  may  note  that 
literature  today  usually  presents  data  plots  using  the  absolute  value 
of  the  amplitude.  One  must  be  careful  when  confronted  by  these  plots. 
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As  Figures  b-b  and  b-ba  have  shown,  the  absolute  value  of  the  amplitude 
does  not  usually  tell  the  whole  picture.  Figure  b-b  plots  the  absolute 
value  of  the  amplitude  in  Figure  b-4.  Figure  b-5a  is  the  contour  plot 
of  Figure  b-b.  The  major  axis  cross-section  is  no  longer  as  well 
defined  by  the  peaks  as  in  Figure  5-4  or  b-4a.  Nevertheless,  producing 
an  image  using  the  spatial  impulse  response  is  very  promising. 

Although  this  thesis'  imaging  theory  is  not  as  'rigorous'  as  Lewis 
and  Bojarski's  work,  the  theory  is  more  flexible  in  application.  The 
object  is  only  required  to  be  illuminated  at  the  aspect  angles  over  a  2* 
solid  angle.  The  shadowed  side  information  is  also  employed  in  the 
image  reconstruction  process.  There  is  no  need  for  any  assumption  on 
the  object's  shadowed  side  geometry.  If  the  start  and  the  end  of  the 
object's  one  dimensional  impulse  response  for  every  aspect  angle  over 
half  of  the  4n  solid  angle  are  well  defined  by  T-j  ' s  (i  1  luminated)  and 
Tf's  (shadowed)  as  described  ear’.er,  then  an  image  of  the  object  may  be 
produced  using  the  T^'s  and  Tf's.  For  Tj's  and  Tf's  not  defined 
distinctly,  careful  interpretation  on  the  spatial  impulse  response  must 
be  used.  In  the  case  of  a  full  4tt  illumination,  the  image  may  be 
reconstructed  using  only  T^'s  information.  This  converges  to 
Lewi s-Bojarski ' s  identity.  This  theory  is  better  in  utilizing  data 
information  available  but  it  lacks  a  concrete  proof. 

Let's  investigate  this  imaging  possibility  further  using  a  six  inch 
diameter  metallic  sphere  whose  centre  is  located  three  inches  off  the 
centre  of  the  reference  plane  on  the  xi-axis.  The  cubic  sampling 


lattice  used  has  a  safety  factor  of  three.  The  frequency  range  used  is 
0  to  12  Ghz.  Again,  the  frequency  data  are  taken  out  of  the  Mie 
solution.  The  sampling  lattice  is  defined  by  Equation  (4-b)  over  half 
of  the  plane  (w<9<2ir).  The  data  are  weighted  by  the  two  dimensional 
cosine  tapering  function.  The  two  dimensional  version  of  Equation  (4-2) 
is  replacing  the  variable  n  by  the  radial  distance  from  the  centre  of 
the  k-space.  The  shape  of  the  function  is  Figure  4-1  rotated  around  the 
weighting  axis.  The  weighted  data  are  Fourier  transformed  into  the 
spatial  domain  and  presented  as  Figure  b-6.  Figure  b-6a  is  the  contour 
plot  of  Figure  b-6. 

Again  the  initial  speculars  trace  the  illuminated  side  of  the 
sphere  nicely.  The  Tf's  are  not  as  well  defined  as  the  finite  cylinder 
case.  More  interpretation  work  is  required.  The  valley  on  the  shadow 
side  shows  a  curvature.  This  may  be  indicating  the  back  side  of  the 
object  having  a  curvature.  The  radius  of  the  valley's  curvature  is 
about  3 7T  inches,  which  is  the  equivalent  distance  a  wave  would  creep 
before  scattering  back  in  the  transmitting  direction  (Figure  b-7).  The 
radius  of  the  curvature  created  by  the  initial  speculars  is  about  six 
inches,  which  is  the  2R  distance  travelled  by  the  wave.  The  first  R  is 
the  distance  travelled  to  the  target.  The  second  R  is  the  distance 
travelled  by  the  wave  scattering  back  from  the  target.  Using  this 
information,  the  circle  with  the  three  inch  radius  can  be  formed. 


Of  course,  if  one  uses  data  over  the  full  plane  as  Lewi s-Bojarski 1 s 
work,  then  there  is  no  need  for  the  previously  described  type  of  pattern 
recognition  interpretation.  The  initial  speculars  usually  define  the 
perimeter  of  the  object,  if  it  is  a  smooth  convex  body.  The  perimeter 
indicated  is  only  one  cross-section  of  the  target  on  the  major  axis.  If 
more  cross-sectional  information  of  the  target  is  available,  then  an 
image  of  the  whole  object  may  be  produced.  The  potential  on  the  image 
reconstruction  is  high,  but  more  research  work  is  required  in  the  area; 
particularly  in  the  pattern  recognition  area. 


TABLE  5-1 

AN  EXAMPLE:  TIME  ESTIMATION 


Type  of 
enclosu re 

Measurement 

time 

(mins) 

Interpolation 

time 

(mins) 

Total 

time 

(mins) 

ci  rcular 

5.42 

40.09 

45.51 

rectangular 

2.32 

8.47 

10.79 

squared 


5.80 


5.80 


rniain.  Data  are  taken  over  the  half 


method 


5-2a.  Contour  plot  of  FI 


transformed  discretely  into  the  spatial  doamin.  Data  are 
taken  over  the  half  plane:  n  <  0  <  2tt 


CHAPTER  VI 


CONCLUSIONS 

If  a  signal  is  limited  in  the  spatial  domain  (or  wave  number 
domain),  this  signal  is  sufficiently  characterized  by  its  values  over 
the  discrete  sampling  lattice  in  the  wave  number  domain  (or  the  spatial 
domain).  In  the  two  dimensional  case,  the  sampling  locations  in  the 
wave  number  space  are  specified  by  the  vector: 

1 1 1  u  i  +  1 2 

where 

llt  12  =  0,  ±  I,  ±  2,  ±  3t  ... 

The  ux  and  u2  are  related  to  v^  and  v2  by  the  following: 

CUi|u2]  =  2Tr[Vl|v2]  T 

where 

-T  :the  transpose  of  the  inverse  of  the  matrix  formed 


The  v^  and  v £  are  specified  by  one's  choice  on  how  the  space  limited 
signal  is  assumed  to  repeat  itself  in  its  domain,  is  pointed  to  the 
centre  of  one  of  the  closest  images  in  the  periodic  lattice.  Vg  is 
independent  of  v^,  and  it  points  at  the  centre  of  another  close  image. 
By  defining  the  settling  time  as  the  end  of  the  impulse  response,  one 
has  a  space  limited  three  dimensional  impulse  response  signal.  In  the 
finite  circular  cylinder  data  presented  in  Chapter  V,  the  signal  is 
assumed  to  be  4  times  the  size  of  the  cylinder  in  the  two  dimensional 
spatial  plane.  The  choice  of  4  is  equivalent  to  choosing  the  one 
dimensional  impulse  response  of  having  a  settling  time  four  times  as 
long  as  the  wave  would  travel  over  the  length  of  the  object  at  any 
aspect  angle.  Or  the  safety  factor  is  chosen  to  be  one.  (i.e., 

2(1+K)  =  4  <=>  K  =  1)  This  factor  provides  sufficient  results  in  both 
Chapter  IV  and  V. 

In  Chapter  V,  the  example  of  the  finite  circular  cylinder  with 
dimensions:  6"  in  length  and  3"  in  diameter,  the  two  dimensional 
containment  unit  is  chosen  to  be  a  square.  The  squared  repetitive 
lattice  in  the  spatial  domain  is  defined  by, 

*1  ■ « is] 

where 

R  =  6" 

and  the  safety  factor  K  is  chosen  to  be  1.  With  these  inputs  to  the 
program  PTGRIO  (see  Appendix  B),  the  program  outputs  the  sampling 
lattice  defined  by  u1  and  uo.  The  data  output  is  arranged  with  the 


aspect  angles  and  the  corresponding  frequency  increment.  One  makes  the 
necessary  data  extraction  at  the  proper  aspect  angles,  and  increments 
the  frequency  until  the  highest  frequency  is  reached.  One  could  also 
have  defined  a  rectangular  repetitive  lattice.  Then 

V1  =  R  U  V1  =  R  [0<5] 


where 


R  =  6" 

Or,  for  a  circular  repetitive  lattice, 

where 

/  2  2 
R  =  /(6")  +  (3" ) 


=  6.708" 

Again  these  Vj^  and  v2  values  can  be  input  into  the  program  PTGRIO  (see 
Appendix  B)  to  obtain  the  sampling  lattice. 

Thus,  one  can  sample  discretely  and  interpolate  in  the  wave  number 


space  to  reproduce  the  frequency  response  of  an  object.  Fourier 
transforming  this  wave  number  signal  into  the  spatial  domain  gives  a 
representation  of  the  spatial  impulse  response  of  the  object.  The 


spatial  impulse  response  is  defined  as  the  image  function  obtained  by 
three  dimensional  Fourier  transforming  the  far  field  frequency  response 
of  a  finite  object  at  all  aspect  angles  defined  over  the  4*  solid  angle. 
Since  most  objects  have  different  shapes,  sizes,  and  orientations,  their 
corresponding  spatial  impulse  responses  are  different.  Their 
corresponding  sampling  lattices  will  also  be  different.  One  general 
sampling  lattice  that  is  applicable  to  a  set  of  objects  is  desirable. 
This  is  accomplished  by  introducing  different  types  of  canonical 
containment  units  to  confine  the  spatial  impulse  responses.  Two  common 
canonical  cells  are  parallelepiped  and  sphere.  Two  dimensional  examples 
are  shown  in  Figures  2-b  and  2-7,  with  their  corresponding  sampling 
lattices  shown  in  Figures  2-6  and  2-8  respectively. 

In  Chapter  IV,  a  six  inch  metallic  sphere  is  chosen  as  an  example 
to  compare  two  types  of  sampling  lattices  -cubic  and  isotropic.  The 
comparison  is  performed  on  the  interpolated  one  dimensional  impulse 
responses  at  different  aspect  angles  using  the  interpolation  scheme 
defined  in  the  sampling  theorem.  As  one  expects,  the  results  turn  out 
to  be  competitive  for  the  two  types  of  sampling.  Efficiencies,  in  the 
sense  of  the  least  number  of  sampling  points,  are  different  for 
different  canonical  containment  cells  used  on  the  same  object. 

Efficiency  in  using  one  type  of  canonical  containment  cell: 

n  =  CONSERVATIVE  ESTIMATE  OF  THE  OBJECT'S  VOLUME 
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Efficiencies  among  a  sphere  cap  cylinder,  a  sphere  and  a  cube  using 
cubic,  isotropic,  and  rectangular  box  confinement  units  are  presented  in 
Chapter  IV.  The  efficiency  definition  proves  to  be  a  very  good  concept 
in  deciding  the  type  of  sampling  lattice  for  an  object  or  a  group  of 
objects.  This  is  under  the  assumption  that  the  settling  time  of  the 
spatial  impulse  response  is  shaped  similarly  to  the  object;  e.g.,  the 
settling  time  of  the  impulse  response  of  a  sphere  is  the  same  in  every 
aspect  angle. 

To  obtain  an  approximation  to  the  spatial  impulse  response  from  the 
sampled  data  over  a  finite  frequency  range,  one  can  use  the  discrete 
Fourier  transform.  Because  of  today's  digital  computer  design,  the 
different  sufficient  characterization  cannot  be  readily  processed 
without  interpolation;  except,  of  course,  the  cubic  lattice  data  sets. 
The  interpolation  step  which  most  people  like  to  avoid,  is  very  time 
consuming.  This  may  be  referred  back  to  the  time  estimation  example  on 
a  sphere  cap  cylinder  presented  in  Chapter  V.  The  interpolation  step  is 
another  factor  that  affects  an  engineering  decision.  The  avoidance 
helps  Mensa  et  al .  to  arrive  at  the  time  response  faster.  The  price 
they  paid  is  the  limitation  of  their  method's  application  to  two 
dimensional  Fourier  transform.  Their  approach  is  thus  not  recommended 
because  of  its  inability  to  be  expanded  into  higher  dimensions.  The 
sampling  criteria  accompanying  their  method  are, 
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the  angular  increment: 


(  <  2D 


X  «  2D 


.  _i  X 

A9  <  <  sin  (^CT)  X  <  2D 


IT 

v  <  T 


2D 


the  frequency  increment: 
c 

A  f  <  2(l+K)D 

where 

D  =  maximum  dimension  of  the  object 
X  =  wavelength  of  the  frequency  used 
c  =  speed  of  light 
K  *  some  safety  factor 

In  the  finite  circular  cylinder  example, 

D  =  17.04  cm 
K  =  0.7b 
A  f  <  0.b03  Ghz 

at  the  highest  frequency  of  10  Ghz, 

X  =  3  cm 
A  9  <  5  degrees 

Therefore, 

the  frequency  increment  chosen  =  O.b  Ghz 
the  angular  increment  chosen  =  1  degree 


Nevertheless,  processed  results  from  UTD  solution  on  a  finite  circular 
cylinder  support  the  convergence  of  Mensa  et  al.'s  method  to  tvo 
dimensional  discrete  Fourier  transform. 

Lewi s-Bojarski 1 s  identity  requires  either  the  object  be  illuminated 
at  all  angles  or  assumption  be  made  on  the  shadowed  side.  Results  on  a 
finite  metallic  circular  cylinder  and  a  metallic  sphere  indicate  that 
the  spatial  impulse  response  approach  does  not  have  the  above 
restriction,  though  proper  interpretation  may  be  required.  Furthermore, 
the  use  of  the  spatial  impulse  response  to  imaging  can  converge  to 
Lewi s-Bojarski ' s  results.  There  is  also  an  indication  that  the 
presentation  in  the  form  of  the  absolute  value  of  the  amplitude  does  not 
necessarily  provide  the  proper  picture  for  pattern  recognition.  The 
plain  amplitude  representation  with  positive  and  negative  values  is 
sometimes  more  appropriate.  This  is  concluded  by  comparing  the 
Figure  b-4,  or  5-4a  with  5-5  or  5-5a.  The  perimeter  of  a  major  axis 
cross-section  of  a  finite  circular  cylinder  is  shown  more  distinctly 
using  the  plain  amplitude  presentation.  Judging  from  the  two 
dimensional  impulse  responses,  one  can  deduce  the  substantial  potential 
of  the  spatial  impulse  response  in  image  reconstruction. 

The  spatial  impulse  response  has  numerous  applications  including 
target  identification  and  imaging.  Although  smooth  convex  metallic  body 
examples  are  considered  here,  tomographic  applications  on  other  types  of 
bodies  are  possible.  In  all,  the  N  dimensional  sampling  theorem 
provides  new  insights  into  the  sampling  criteria  in  the  wave  number 
space  for  a  finite  object.  The  potential  in  reduced  management  time  on 


CHAPTER  VII 


RECOMMENDATIONS 


Traditionally,  two  and  three  dimensional  data  are  presented  in 
cubic  lattices,  for  which  present  day  digital  computers  are  designed. 
Although  the  digital  computer  of  today  manages  data  in  cubic  lattices, 
the  cubic  lattices  are  not  necessarily  the  most  efficient  in  terms  of 
the  least  number  of  data  samples.  The  most  general  approach  to  solve 
this  computer  problem  requires  the  cubic  data  management  structure  of 
the  digital  computer  be  modified  into  a  more  general  data  structure.  In 
another  words,  data  organized  in  any  random  fashion  can  be  processed  by 
this  computer.  If  this  general  approach  is  not  practical,  the  next  best 
step  is  a  faster  interpolation  scheme  either  in  hardware  or  software. 

The  lowest  level  on  the  hierarchy  of  improvement  is  the  improvement  for 
specific  application.  In  this  thesis,  a  general  Fourier  transform  that 
can  perform  on  any  data  lattice,  fits  into  this  category. 

After  some  of  the  computer  problems  are  solved,  the  next  step  in 
the  development  is  to  account  for  the  experimental  noise.  How  does  one 
extract  the  true  information  that  is  embedded  in  noise?  Without  this 


generalization,  this  report  is  only  useful  in  information  storage  and 
retrival  on  noise  free  data.  Noise  free  is  in  the  sense  that  the  noise 
effect  in  measurement  is  eliminated  before  storage.  The  full 
development  of  this  report's  theory  will  enable  significant  reduction  in 
time  on  data  measurement,  processing  and  storage. 

The  whole  report  is  focused  on  the  monostatic  data.  It  would  be 
interesting  to  see  how  this  theory  holds  up  with  the  bi static  data. 

Using  the  definition  of  the  three  dimensional  Fourier  transform,  one  can 
derive  a  similar  method  that  parallels  Mensa  et  al.'s  approach. 

By  converting  the  k-space  coordinates  into  spherical  coordinates: 

kj  =  psin9coS(j>  k ^  -  psin9sin<j>  k^  =  pcos0  (7-1) 

equation  (1-1)  becomes: 

1  -  Jlxjpcos  <  (x,  p)  2 

f(x)  =  I  I  J  Hp»  e»  4>)e  p  sin0d<f>d9dp 

p=0  0=0  <j.=0  (7-2) 


where 

M  =  /cf  +  x2  +  x§ 

(Rotate  the  k-space  coordinate  system  so  that  the  kj-axis  is  in  line 
with  x  :  =>  <  (x,  p)  =  0). 


Then, 


=  8 7S  $  I 


-  j  |  x  |  pco  s  9  o 

/  F(p,  6,  <t»)e  P  Sin9d4>d9dp 


p=0  9=0 


(7-3) 


By  considering  one  particular  aspect  angle:  4>m,  9^ 


sin9^  ^  -  j  |  x |  pcos £ 

fimU)  =  "sTT  J  F(P.  9i.  $m)e  p  dp 


(7-4) 


With  all  aspect  angular  data. 


f(x-)  -  I  l  f1m(x-) 
i  =1  m=l 


(7-t>) 


=  the  impulse  response  of  the  finite  object  at  x 


Even  though  Equation  (7-4)  requires  the  integration  over  all 
frequencies,  the  integral  can  be  approximated  over  three  regions:  the 
Rayleigh,  the  resonance,  and  the  optical.  Thus,  the  integral  in 
Equation  (7-4)  is  not  impossible  to  be  solved.  There  are  problems 
associated  with  this  approach.  The  Nyquist  angular  requirement  is 
frequency  dependent  (Equation  (2-17)).  In  order  to  satisfy  the  Nyquist 
angular  requirement  at  high  frequency,  the  signal  is  excessively  sampled 
at  the  low  frequency  spectrum,  but  sampled  only  adequately  at  the  high 
frequency  spectrum  for  a  finite  frequency  range  of  interest. 

Nonetheless,  the  approach  is  viable  if  one  does  not  intend  to  extend  the 
bandwidth  of  the  approximated  spatial  impulse  response. 
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With  the  help  of  the  different  sampling  criteria,  the  infinite 
number  of  samples  required  to  reconstruct  the  impulse  response  over  a 
finite  frequency  range  has  changed  to  a  finite  number.  Though  the 
number  is  finite,  the  measurement  time  can  be  very  substantial  when  an 
expanded  frequency  range  in  two  dimensions,  or  three  dimensions  is 
required.  Professor  Leon  Peters  suggested  another  approach  to  further 
reduce  the  measurement  time  [15],  At  high  frequency,  a  target's 
frequency  response  is  mostly  contributed  by  its  major  scattering 
centres.  If  one  can  make  a  set  of  different  canonical  scattering  centre 
measurement,  then  most  targets'  high  frequency  response  can  be  built 
using  the  proper  phase  shift  factors.  For  most  scattering  centres, 
their  frequency  responses  are  relatively  simple.  As  a  result,  the 
Nyquist  criteria  for  these  centres  in  the  frequency  domain  are  more 
relaxed  than  complete  structures.  These  canonical  scattering  centre 
data  can  be  reused  to  reconstruct  the  high  frequency  response  of  other 
targets.  In  another  words,  after  the  canonical  scattering  centre  data 
are  available,  one  only  measures  the  low  frequency  spectrum  before  the 
one  dimensional,  two  dimensional,  or  three  dimensional  impulse  response 
of  an  object  can  be  reconstructed.  This  would  be  another  interesting 
area  for  further  exploration.  However,  more  development  work  is 
required  in  all  these  described  areas  to  extend  this  report  into  a  more 
useful  engineering  tool. 
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APPENDIX  A 


THE  DERIVATION  OF  EQUATION  (2-13) 


( 

2“l  cos  \7T" )  cos  (Rw2) 


G  (  wi  »  “o )  =  2  2  2  X  \ 

1  1  R\(u>J  -  3002)  ' 


f  2RwA 
-2^  cos  \7J~ ) 

/r»A 

-2/T  u2  sin\7Ty  sin  (Ro>2 


Equation  (65)  of  Petersen  and  Middleton  [7]: 


Regular 

Hexagon 


Equation  description  for  the  regular  hexagon:  (see  Figure  A-l) 


© 


X2  =  -/3  xi  +  2R 


gure  A-l.  Hexagonal  integration  limit 
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'  x2  =  /3  XI  -  2R 

f  r\  ( 2R\ 

\rsj  <  xi  <  \73y 


G(-)  =  — L_  ^  +  l2  +  l3) 

2/3 R 


where 
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ii  ■  / 


j  (“lxi  +  0)2x2) 

/  e  dx2  dxi 


«!  .(H)  ^  -® 


© 


1 2  -  / 


xi 


■00 


j (“lxi  +  U)2  X2) 

/  e  dx2  dxi 

x2  =-R 


1 


(A-l ) 


(A-2) 


(A-3) 
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j“2 


< 


j(“l  +  /3o>2)xi  +  2jRa»2 


j(“l  +  /Tui2) 


/3 


-2R 

71 


Therefore, 


j(wl  -  /3w2)xi  -  2jR<*>2 


j{“l  -  /3ui2) 


-R 

/3 


-2R 

71 


2jRa>2 


*1  " 


u>2(<*>1  +  /3«>2) 


-j(a)i  +  /3u>2)(2«)  -j(a>!  +  /3u)2)(i.) 

73  71 


-  e 


-2jRa)2 


u2(wi  -  /3to2) 


-j(a)i  -  *'3(D2)(^.)  -  /3a>2)(2K) 


-  e 


Equation  (A-3): 
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R 

,  j(“lxi  +  u)2x2) 

J  e 


-R  x2  =  -R 
X1  =  /3 


R 


>1X1 

e 

rs 

>2*2 

e 

>1 

-R 

>2 

rs 


(J\  Ai«  'N 

\»1 J  si n  \V3 J 


sin  (u>2R ) 


Therefore, 


Equation  (A-4), 


13  -  I 


-/lx 1  +  2R 

/ 


j( wixi  +  “2X2) 

e  dx2dxi 


R 

xi  =1/3 


X2  =  ^3xi  -  2R 


j“lxi 


X1 


j“2X2 


J“2 


-/lx 1  +  2R 


dxi 


/3xi  -  2R 


j“lxi  f  jo)2(-/3xi  +  2R)  jo>2(/3xi  -  2R) 


J“2 


-e 


dxi 


xi  =1/3 


2R 
71, 

_L  \  f  j(wl 

j“2  *  [e 

*1  ih) 


/3o>2)xi  +  2jRo>2  j  ( o>i  +  /3o>2)xi 
-e 


2jRo>2 


dxi 
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Ju2 


j(wl  -  /3o)2)xi  +  2jKu>2 


j(“l  -  /Ti>2) 


R 
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j(wl  +  73<u2)xi  -  2jRo)2 
e 

71 
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R 

J? 

Therefore, 
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-  /3u2)  L 
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j(«l  +  /3u2)(2k)  j(wl  +  *'3u2)(^.) 
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2cos 


(“1  +  /3u)2)  R  -  2Ro>2 
77 


(This  expression  continues  on  the  next  page.) 


2Ruj2] 
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2Ru)2J 
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[Since:  cos  <*  -  cos  6  =  -  2  sin  7  (a  +  6)  sin  J  (a  -  8)] 
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rs  ~r 


“2(^1  +  /3u)2) 


“ATJ 
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Equation  (A-l), 


G(^)  = - 2  +  *2  +  13) 

2/3R 


r  o>i  ( a>i 


'  /  Ru)l\ 

2  <*>i  cos(  “73  J  cos  ( Ro>2 J 

/2R<*>]\ 

-  2  oil  cos  (  73  J 

/RioiN 

-  2  /I  o>2  sin^~7I J  sin  (R<*>2) 


=  Equation  (2-13). 
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Figure  B-l.  Flow  chart  showing  the  relationship  among  the  written 
programs 
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PROGRAM; 


RECCNSHD 


THIS  PROGRAM  WILL  TAKE  A  DATA  SET:  'CDEFF'  IN  REAL  AND  IMG  INARY  FORMAT 
AND  INTERPOLATE  THE  POINTS  IN  BETWEEN  SAMPLES  USING  SAMPLING  THEOREM 
APPROACH.  PRESENTLY,  THIS 

PROGRAM  IS  ALLOWED  TD  TAKE  100  SAMPLES.  IF  FORE  IS  REQUIRES,  PLEASE 
CHANGE  THE  ARRAY  SIZE  OF  CDEFF  AND  THE  VALUE  OF  NUM. 

MSIZE:THE  ARRAY  SIZE=2*  (NUFCEH  OF  SAMPLES)  *1 
NUMsNUMBER  OF  SAMPLES 
XINC:  INCREMENT 

THE  SAMPLES  ARE  TAKEN  AT  DELTA  FREQUENCY  OF  'FRED*. 

THE  EREDUENCY  RANGE  FOR  INTERPOLATION  IS  SPECIFIED  BY  '  INITIAL' 

AND  'LAST'  -  THE  NUM3ER  OF  POINTS  IN  'DATA*  IS  'NPl*.  THE 
OUTPUT  FILE  IS  SPECIFIED  BY  'FNAME2'. 

THIS  PROGRAM  WILL  LINK  WITH  SINC,  JPLOT,  SUM1D  AND  'PLOTLIB 

REAL  INITIAL.  LAST 
COMPLEX  CDEFF (201)  .SUM, DATA (200) 

CHARACTER*10  FNAMEl -FNAEE2 
CHARACTER  0DM1 
NUM100 
MSIZE>2*NUM  +1 
NSTEP^200 
DHTIAL-6.E7 
LAST*12.E9 

XINO  (LAST-INITIAL) /(NSIEP-1) 

PI-3.14159265 

THE  SUFtMTON  IS  FROM  -N  TO  +N 

CDEFF  -  THE  ARRAY  OF  COEFFICIENTS  FOR  iTH  TERM  IN  THE  SUFMCTTON 
WHERE  THE  INDEX  -N  IS  RES  PRESENTED  BY  1 

+N  IS  RES  PRESENTED  BY  2*N  +1 
T  -  ASSUFED  CUTOFF  TIME-  1/ (2*SAMPLING  FREQUENCY) 

DATA  STRUCTURE  OF  FNAFCl 
N  :NUFCER  OF  POINTS  IN  THE  FILE  (•) 

EMIN  :SF«LLEST  EREDUENCY  USED(*) 

FRED  .SAMPLING  FREQUENCY  (*) 

CDEFF  ;i)  s  (REAL  -  IMG  INARY)  (*) 

TYPE  *.  'TYPE  FILE  NAAE  CONTAINING  COMPLEX  CDETFIONT5' 

ACCEPT  2 -FNAMEl 
2  FORMT  (A10) 

OPEN  (UNT>8,NAME> FNAMEl. TYPE- 'CLD'> 

READ (8,*)  N 
READ(8,*)  EMIN 
READ(8.*)  FREQ 
T>1./(2.*FRBD) 

IF  (N.GT.FDM)  GQ  TO  8888 


V-V- 


•-»  »  •  »  *»  •  -  *  »  •  *  v* 

■ .... 


133 


i 

i 


NDIM-2*N  +1 
C 

C  INITIALIZATION 

C 

DO  5  I-1.MSI2E 

CD  EFF(I)  •(£)., 0.) 

5  CONTINUE 

READ (8,*)  CDEFF(NUM) 

DO  10  I-l.N-1 

READ(8,*)  ODEFFUWW-I+l) 

C 

C  FILL  IN  COMPLEX  CONJUGATE  FOR  THE  NEGATIVE  FREQUENCY 

C 

CDEFF(HJMfl-I)*CDNJG(OOEFF(NUM+I+l) ) 

10  CONTINUE 

CLOSE  (UNITES  ,DISP* '  SAVE 1 ) 

IEEROO 

C 

C  SET  ANY  FREDUENCY  INTERPOLATION  SJECIFTED  BELOW  FMIN  TD  ZERO 
C  NZERO:  NUFBER  OF  ZEROS  TD  BE  PLACED 

C 

IF  (ETON.  LE.  FREQ)  GO  ID  127 
FC  ERO=  INT  (ETON/XINC) 

DO  123  I-l.FEFRO 

DATA(I)«(0..0.) 

123  CONTINUE 
C 

C  STARTS  THE  INTERPOLATION 

C  (THE  SUWftTIDN  IS  PERFORMED  BY  FUNCTION  SUM1D) 

C 

127  ILAST"INT  ( (N-l )  *FREQ/XINC) 

IF  (ILAST.LT.  NSTEP)  INDEX-ILAST 
IF  (ILAST.GE.  NSTEP)  INDEX ^STEP 
DO  20  I-fEERO  +1.  INDEX 

DATA  (I)  -  SUMID  ( ( (1-1 )  *XINC+ INITIAL)  *T*2 .  *PI ,  CDEFF.  MSIZE.  NDIM) 
20  CONTINUE 

C 

C  SET  ANY  SPECIFIED  FREQUENCY  INTERPOLATION  HIGHER  THAN  FTONHWRBQ 
C  TD  ZERO 

C 

DO  133  I« INDEX  +1. NSTEP 
DATA(I)-(0.,0.) 

133  CONTINUE 

C 

C  PLOT  AND  OPTIONAL  WRITE  DJID  A  FILE:  FNAFC2 
C 

CALL  J  PLOT  (DATA.  NSTEP,  INITIAL. XINC, 0 .NSTEP-1) 

23  WRITE (6,*)  'WRITE  THE  INTERPOLATED  DATA  INTO  A  FILE?’ 

WRITE  (6,*)  'DYES,  IN  AMPLITUDE  AND  PHASE  (RADIAN)  FORM’ 

WRITE (6,*)  ’2) YES,  IN  REAL  AND  IFW3INARY  FORM' 

WRITE  (6,*)  '3>fO.  FORGET  IT  I ' 


vvv 


ACCEPT  25/00. 

25  FORMAT (Al) 

IF  (COK1.BO. '3M  GO  TO  9999 

IF  ((COMl.tE.  '1')  .AMD.  (00M1.NE.  '2'))  GO  TO  23 

WRITE (6.*)  'STORAGE  FILE  NWC: ' 

ACCEPT  30,FNAME2 
30  FORMAT  (AlO) 

IF  (0»0. ED. '2')  00  TO  35 
C 

C  CONVERT  TO  AMPLITUDE  (LINEAR)  AND  PHASE  (RADIAN) 

C 

DO  33  I-l.NSTEP 

XREAL-CABS  (DATA  (I ) ) 

XIMAG-AIAN2  (AIMAG  (DATA (I ) ) , REAL  (DATA (I ) ) > 

DATA(I)  -OULX  OCREAL.XIMAG) 

33  OONTENUE 

C 

C  OUTPUT  TO  FNAW2 

C  STRUCTURE  OF  PNAtC2: 

C  DATA (I) -.FREE  (DMPLEX  FORMAT (*) 

C 

35  OPEN  (UNIT-8. NAfE«FNAME2.TXPE» 'NEW') 

DO  40  I-l.NSTEP 

WRITE (8.*)  DATA (I) 

40  CONTINUE 

ODSE  (UNIT-8,  DISI^' SAVE') 

GO  TO  9999 

8888  WRITE (6.*)  'ERROR: NOT  ENOUGH  SPACE  SPECIFIED  IN  THE  PROGRAM' 
9999  STOP 

END 
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SUBROUTINE 


JHOT 


THIS  SUBROUTINE  WILL  DO  RECTANGULAR  PLOT  ON  A  COMPLEX  ARRAY:  'DATA' 
FDR  REAL  AND  IMAGINARY  PLOT  OR  EftGNITUDE  (LINEAR)  AND  PHASE  (RADIAN) 
PLOT. 

DATA  :  ARRAY  NAFC 

NPOINT  :THE  ARRAY  SIZE 
EMIN  .-SMALLEST  ELEMENT  OF  THE  ABSdSSDR 
FIN®  :THE  INCRMENT  SIZE 
NSTART  :TH£  START  PLOTTING  INDEX 
NLAST  :THE  LAST  PLOTTING  INDEX 
IF  (MAST. GT.  NPOINT)  THE  LAST  POINT  PLOTTED  IS  NPOINT 

THIS  REQUIRES  THE  SUPPORT  OF  'ILOTLIB. 

SUBROUTINE  JPLOT  (DATA,  NPOINT,  EMIN .  FIN®.  NSTART,  MAST) 

CHARACTER  OOM 
COMPLEX  DATA  (NPOINT) 

DIMENSION  YAXISl (9000) ,YAXIS2 (9000) ,XAXIS (9000) 

PI  >3. 14159265 
MSIZB=9000 

INITIAL  IZATDON 

DO  3  I-l.MSIZE 

YAXISl  (I)  *0. 

YAXIS2(I)»0. 

XAXIS(I)-0. 

CONTINUE 

WRITE (6.*)  'DO  YOU  WANT  REAL  AND  IMAGINARY  PL0T7Y/N' 

ACCEPT  6,  COM 
FORMAT  (Al) 

IF  (NLAST. GT. NPOINT)  MAST»NP0INT 
NP-NLAST-NSTART+1 
IF  (OOM.  ED.  'N')  GO  TD  15 
IF  (OOM.NE.  'Y' )  GO  TO  5 

REAL  AND  IMY3INARY  PREPARATION 

DO  10  1*1  .NP 

YAXISl  (I)  -REAL  (DATA (I+NSTART) ) 

YAXIS2 (I )“AIFWG (DATA (I+NSTART) ) 

XAXIS(I)-(I-1+NSTART)*FIN®  +  FMIN 
)  CONTINUE 
GO  TD  25 

GENERATE  AMPLITUDE  AND  PHASE 

15  DO  20  1*1 .NP 

XAXIS  (I)  “(I+NSTART-1)  *FIN®  +FMIN 
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OOMBLEX  JUNCTION 


SUKLD 


C 
C 
C 

C  THIS  SUBROUTINE  WILL  ID  SUEWmDN 

C  X  -  NEW  PARAMETER 

C  Y  -  COEFFICIENT  ARRAY 

C  MS IZ E-SIZE  OF  THE  COEFFICIENT  ARRAY 

C  N  -  NUfBER  OF  COEFFICIENT  ARRAY  ELEMENTS  THAT  HAS  t©N-ZERO  VALUES 

C  H5TE:BOTHE  MSIZE  AND  N  ARE  ODD  NUEBERS 

C 

CDMH-EX  FUNCTION  SUM1D(X,Y,MSIZE,N) 

EXTERNAL  SINC 
COMPLEX  Y  (MSIZE) 

PI*3. 14159265 
MEI>=MSIZE/2  +1 
SUMID-Y  (MED)  *SINC  (X) 

DO  5  I-l.N/2 

SUMlO-SUMlIHY  (MSIZ  E/2-1+1 )  *SINC  0(+I*PI)  +Y  (MSIZ  E/2+1+1)  *SINC(X-I*PI) 
5  CONTINUE 

RETUWJ 
END 
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PROGRAM  : 


PTGRID 


IMIS  PROGRAM  WILL  GENERATE  ALL  THE  ANGLES  AND  FREDUENCLES 
FOR  THE  GRID  IN  THE  FREQUENCY  DOMAIN.  THE  DATA  ARE 
ARRANGED  IN  ORDER  FROM  SMALLEST  ANGLE  TD  THE  LARGEST  ANGLE. 

(0  TD  2*PI)  THE  SMALLEST  FREQUENCY 

RADIUS  IS  STORED.  SO  AT  THE  TIFE  OF  MEASUREMENT,  ONLY 
MJLTIPLICATTON  OF  THE  RADIUS  IS  NECESSARY. 

THE  GRID  iOINTS  ARE  WRITE  INTO  FILE:  'OUT.  DAT' 

THE  DATA  FORMAT  OF  'OUT. DAT': 

FMIN-FMAX  : FREQUENCY  RANGE  SPECIFIED(2E15.8) 

FDELTA  : SAMPLING  FREQUENCY  (E15 . 8) 

NFOINT  :NUFBER  OF  POINTS  IN  THE  FILE(I8) 

ISO, RADIUS  :ISO»N.  LON-ISOROPIC  SAMPLING  IS  USEDCA2) 

:ISO-T,  ISOTROPIC  SAMPLING  IS  USED 

:  RADIUS:RADIUS  OF  THE  ISOTROPIC  CELL  IN  WKE15.8) 

V1.V2  : VECTORS  USED  TD  DEFINE  THE  SAMPLING  LATTICE (4E1 5. 8) 

U1.U2  : VECTORS  USED  TO  DEFINE  THE  PERIODIC  LATTICE (4E1 5. 8) 

DATA(*,1)  ,DATA(*,2)  :  DATA(4E15.8) 

DATA(*,1)  IS  THE  MEASUREMENT  ANGLE  AND  FREQUENCY  ARRAY 
DATA ( * . 2 )  IS  THE  INDEX  SPECIFYING  HOW  MANY  UNITS  OF 
VI .V2  ARE  USED 

WARNING:  SIZE  OF  DATA  IS  (10000,2) 

LINK  PTC RTO.GEN1.GFN2. SEARCH. INSERT, POSH. NGRMAL.'SSP 

COMPLEX  VI .V2 -V3 .V4 -U1  .U2. DATA (10000. 2) 

DIMENSION  WORKl(2)  ,WORK2(2) 

REAL  MAGI  .MMG2  -MAXI, MARGIN. MATRIX (4) 

CHARACTER  ISO 
MSIZ  E>=10000 
PI«3. 14159265 

START  WORKING 

WRITE  (6.*)  '  LOWEST  AND  HIGHEST  OPERATING  FREQUENCY  IN  HERTZ ' 

ACCEPT  *,  EMDJ.EMAX 
IF  (FMIN.GE.  EWAX)  GO  TD  2 

WRITE  (6.*)  '  DO  YOU  WANT  TD  DO  ISOTROPIC  SAMPLING  ?  T/F  ' 

ACCEPT  4. ISO 
FORMAT  (All 

IF  (ISO. ED.  'F')  00  TD  5 
IF  (ISO. NE. 'T')  00  TD  3 

WRITE  (6.*)  '  DIAMETER  OF  NORMALIZATION  IN  FH?' 

ACCEPT  *, RADIUS 


DEFTNITIDN  OF  ISOTROPIC  VECTORS 


Xl«(90KT(3.))/2. 

Yl— 0.5 
X2-0. 

Y2-1. 

GO  ID  € 

CASE  OF  PARALLSLPIPEDIC  CONFINEMENT 

WRITE  (6.*)  'THE  OBJECT  IS  ASSUfED  TO  BE  CONFINED  TO  A  PARALLELOGRAM  1 
WRITE  (6.*)  '  NORMALIZATION  FACTOR  IN  *M' 

ACCEPT  ‘.RADIUS 

WRITE  (6.*)  '  VECTDR1  IN  X.Y  s' 

ACCEPT  *,  Xl.Yl 

WRITE  (6.*)  '  VECIDR2  IN  X,Y  s’ 

ACCEPT  *  X2-Y2 

WRITE  (6!*)  '  SAFETY  MARGIN  FACTOR  CN  THE  TIME  WAVEFORM' 

ACCEPT  ‘.MARGIN 
UlKMPLXOCl-YD 
U2KJ4PLXCX2.Y2) 

MATRIX  (1)  -(XI) 

MATRIX  (2)  «<Y1) 

MATRIX  (3)  «(X2) 

MATRIX  (4)  -(Y2) 

MINV  WILL  DO  MATRIX  INVERSION.  THIS  SUBROUTINE  IS  IN  SSP 

CALL  MINV  (MATRIX.  2  .D.  WORK1  .WORK2) 

KXTE  THAT  VI  AND  V2  IS  TAKEN  FROM  THE  TRANSPOSE 

OF  THE  INVERSE  OF  MATRIX.  THIS  IS  TD  BE  CONSISTENCE  WITH 

DOT  PRODUCT  OF  U(I)  AND  V(J)  -  DELTA (U) 

U  IS  IN  SPATIAL  DOMAIN 

V  IS  IN  EREOUENCY  DOMAIN 

ALSO  V1»(X-COMPONHJT»Y-CDMFONIOT) 

VI -CMPLX  (MATRIX  (1)  .MATRIX  (3) ) 

V2<MFLX  (MATRIX  (2)  .MATRIX  (4) ) 

MAG1<ABS(VD 

PHASE1 WCTAN2  (ADftG  (VI ) ,  REAL  (VI ) ) 

MAG2«CABS(V2) 

PHASE2-ATAN2 (AIMAG  (V2) .REAL (V2) ) 

THETA-ABS  (PHASE1-IHASE2) 

PHASES  ARE  NORMALIZED  TO  THE  RANGE  0:2*PI 

CALL  NDRWiL  (PHASED 
CALL  NORMAL (PHASE2) 


S FACTOR:  EREQUENY  SAMPLING  FACTOR 


non  non  nnn  oooom  t-  non 


S  FACTOR®  ( 3 .  Ell )  /  (RADIUS*MARG  IN) 

IF  <(SFACTOR.LE.FMIN)  .OR.  (FMIN.EC.O.) )  FDEXTA=SFACTOR 
IF  ( (SFACTOR.  GT.  FMIN)  .AND.  (FMIN.  JE.  0 . ) )  FOELTAfFMIN 
RADIUS®!/  (FDELTA*2  - ) 

MAXI®  (FMAX/FDELTA) 

STORE  THE  FIRST  TWO  VECTORS  INTO  THE  ARRAY  'DATA' 

IF  (PHASE2.LT. PHASED  00  TO  11 
DATA  (1 ,1 )  *<MPLX  (JHASE1  .MN31 ) 

DATA  (1,2)  KMHJC  (1.0) 

DATA(2.1XMPLX(FHASE2.»G2) 

DATA(2.2)  -CMHjX(O.I) 

00  TO  12 

L  DATA(l.l)<MHJ((FHASE2.fBC2) 
nATAQ.2)<MPLX(0.1) 

DATA  ( 2 . 1 XMHJC  (PHASEl .  MA31 ) 

DATA(2-2)=<MFLX<1.0) 

2  NP0INT-2 

TO  GENERATE  THE  FIRST  QUADRANT  POINTS 
POINTS  ARE  GENERATED  AND  STORED  IN  SUBROUTINE  GEN1  AND  GEN2 

IF  (THETA.  LE.P1/2.) 

k  CALL  GEN1  (VI  ,V2 -MAXI. MAXI. DATA,MSIZE, NFOINT.l  .1) 

IF  (THETA. GT. PI/2® ) 

k  CALL  GEN2(V1.V2. MAXI, MAXI. DATA#MSIZE,NPODrr,l.l) 

IF  (NPOINT.  GT.  MSIZE)  00  TO  688 
V3— VI 

PHASE3  -  ATAN2  (ADBG  (V3 )  ,REAL(V3) ) 

CALL  NORMAL (PHASE3) 

V3KMPLX(PHASE3.CABS(V3) ) 

SEARCH  FOR  LOCATION  :LOC  TO  PLACE  THE  POINT 

LOOSEARCH  (V3  .DATA.  MSIZE.  NPOINT) 

INSERT  THE  POINT  INTO  THE  TOPER  LOCATION 

CALL  INSERT  (DC,  NFO INT,  DATA,  MSIZE.  V3.-l.,0.) 

TO  GENERATE  THE  2ND  QUADRANT  POINTS 

IF  (THETA.  LE.  PI/2. ) 

k  CALL  GEN2(V1 ,V2 -MAXI. MAXI. DATA, MSIZE, NPOINT, -1,1) 

IF  (THETA,  OT,  PI /2.) 

k  CALL  GENl (VI ,V2. MAXI. MAXI. DATA, MSIZE, NFOINTf-l,l) 

IF  (NPOINT. OT. MSIZE)  GO  TO  888 
V4— V2 

FHASE4  -  ATAN2  (AIMAG  (V4)  .REAL (V4) ) 

CALL  NORMAL (EHASE4) 


*.  v  *_  *.  \ 
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V4<MFLX  (PHASE4  .CABS  (V4) ) 

LOO  SEARCH  ( V4  .DATA.MSIZE,  NTODJT) 

CALL  INSERT  (IOC,NFDIOT,DATA,MSIZE,V4 .0 . ,-l.) 

C 

C  NOW  TD  GENERATE  THE  OTHER  TWO  QUAERAOTS 
C 

IF  (THETA.  I£.  PI/2.) 

6  CALL  GEN2(VX.V2.fftXl.«ua.EATA.KSI2E,NP0INT,I,-I) 

IF  (THETr.OT.PI/2.) 

k  CALL  GEM  (V1.V2. MAXI. MAXl»DATA,MSIZE,NFDINT»l.-l) 

IF  (NFOINT.GT.MSIZE)  GO  TD  888 
IF  dHTQLLEtPI/2*) 

k  *  CALL  GEM  (V1.V2  .MAXI, MAXI. OATA,MSIZE,NFOINT,-l,-l) 

IF  (THETA. GT. PI/2.) 

k  CALL  GEN2(Vl,V2.WOa.(AXl.DATA,MSIZE.NFOINT.-l.-l) 

IF  (NFOINT.Gr.MSIZE)  GO  TD  888 
C 

C  OUTHJT  GRID  POINTS  INTO  FILE:  'OUT. DAT' 

C 

OPEMUNir-a.NAME-'aUT'  .TXPE-'NEW') 

WRITE (8.301)  FMIN.FMAX 
301  FORMAT (2E15 .8) 

WRITE  (8,305)  FDELTA 
305  FORMAT  (E15. 8) 

WRITE  (8,310)  NBOINT 
310  FORMAT  (18) 

WRITE  (8,315) ,  ISO, RADIUS 
315  FORMAT  (A2- El  5 -8) 

WRITE (8.320)  VI -V2 
320  FORMAT  (2E15. 8,213. 5. 8) 

WRITE (8,320)  U1.U2 
DO  20  J-l.NFOINT 

WRITE  (8,320)  DATA(J,1)  ,DATA(J,2) 

20  CONTINUE 

CLOSE (UNIT-8, DISF*' SAVE') 

GO  TO  999 

888  WRITE (6.*)  ' ERROR : S lECI FIED  ARRAY  SIZE  TDO  SMALL  1 ' 

999  STOP 


uou  uuu  ouu 


SUBROUTINE 


GENl 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  SUBROUTINE  WILL  GENERATE  ALL  THE  GRID  POINTS  WITHIN  THE  TWO 
VECTORS  VI .V2  ON  THE  TWO  DIMENSIONAL  IUANE.  THAT  IS  IF  THE  ANGLE 
IS  ACUTE. 

VI  : FIRST  VECTOR 
V2  :  SECOND  VECTOR 

MAX  : MAXIMUM  FREQUENCY  IN  EACH  VECTOR  DIRECTION 

AMAX  : ABSOLUTION  MAXIMUM  FREQUENCY 

DFILE  :DATA  ARRAY  FOR  STORAGE 

MSIZE  :SIZE  OF  DFILE 

M  : NUMBER  OF  POINTS  GENERATED 

NISIGN  :SIGN  OF  Nl  (N1*V1) 

N2SIGN  :SIGN  OF  N2  (N2*V2) 

THIS  REQUIRES  THE  SUPPORT  OF  SEARCH- NORMAL. INSERT 
SUBROUTINE  GENl  (VI  ,V2. MAX.  AMAX, DFILE,  MSIZE,M,NLSIGN,N2SIGN) 
EXTERNAL  SEARCH 

OOMELEX  DFILE  (MSIZE.  2)  .VECTOR.  VL  ,V2 
REAL  MAG, MAX 


CALCULATE  THE  NUMBER  OF  UNITS  IN  EACH  VECTOR  DIRECTIONS 


N1«INT(MAX/(CABS(V1)))  +1 
N2»INT  (MAX/ (CABS  (V2) ) )  +1 
DO  5  1-1. N2 

DO  10  J-1.N1 

VECTORS  *V1*N1SIGN  +I*V2*N2SIGN 
MAG-CABS(VECroR) 

IF  WG.GT.AMAX)  GO  TO  5 
IHASE=ATAN2  (AIMAG  (VECTOR)  .REAL  (VECTOR) ) 

CALL  NORMAL  (PHASE) 

VECTDR=®HJC  (PHASE,  MM3) 

SEARCH  THE  LOCATION 

LOOSEARCH  (VECTOR.DFILE, MSIZE, M) 

RJ-FIOAT(J)  v* 

RI»FLQAT(I) 


INSERT  THE  DATA 


t 

& 

& 


IF  ((NlSKN.LT.O).AND.  (N2SIGN.LT.0) ) 

CALL  INSRT  (LOC.M.DETLE, MSIZE, VECTOR. -RJ.-RI) 
IF  ( (NISIGN.  LT.O)  .AND.  (N2SIGN.CT.0)  > 

CALL  INSERT  (LOC.M, DFILE, MSIZE, VECTDR.-RJ.RI) 
IF  ((N1STON.GT.0)  .AND.  (N2SIGN. LT.O)) 

CALL  INSERT  (IOC,  Mi,  DFILE, MSIZE.VECTOR.RJ.-RI) 
IF  ( (NISIGN. GT.O)  .AND.  (N2SIGN.GT.0) ) 
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C 

C  SUBROUTINE:  GEN2 

C 

C  TBIS  SUBROUTINE  WILL  GENERATE  POINTS  BETWEEN  VECTORS  VI  ,V2  IF 
C  TOE  ANGLE  BETWEEN  THEN  IS  MORE  THAN  90  DEGREES  BUT 
C  LESS  THAN  180  DEGREES 
C 
c 


c 

VI 

:  FIRST  VECTOR 

c 

V2 

:  SECOND  VECTOR 

c 

MAX 

sMAXHUM  FREDUENCY  IN  EACH  VECTOR  DIRECTION 

c 

AMAX 

: ABSOLUTION  MAXIMUM  FREQUENCY 

c 

DFILE 

:DATA  ARRAY  FOR  STORAGE 

c 

MSIZE 

:SIZE  OF  DFILE 

c 

M 

:NUWER  OF  POINTS  GENERATED 

c 

N1SIGN 

:SIGN  OF  Nl  (N1*V1> 

c 

N2SIGN 

:SIGN  OF  N2  (N2*V2) 

c 

C  THIS  REQUIRES  THE  SUPPORT  OF  GEN1 .INSERT, SEARCH 
C 

SUBROUTINE  GEN2(V1.V2.MAX.AM0C.DFILE,MSIZE»M,NLSIGN.N2SIGN) 
EXTERNAL  SEARCH 

CDMH-EX  DFILE  (MSIZE.  2) , VECTOR, VI  ,V2 
REAL  MAG, MAX 
LOGICAL  SAW! 

PI-  3.14159265 

THETft-ABS  (ATAN2  (AIMAG  (VI  > , REAL  (VI ) )  -ATAN2  (AIMAG  (V2)  ,REAL  (V2) ) ) 
XMAX-MUC/SIN  (THETA) 

CALL  GEN1  (VI  .V2  .XMAX,  AMAX. DFILE.  MSIZE, M.NISIGN.EESIGN) 

IF  (THETA. ED. FI/2.)  00  TO  555 
(QCDRRECM) 

N2CDRRECT-0 

IF (AMOD (XMAX, CABS (VI) )  .GT. 0.)  NlOORREOl 
IF (AMOD (XMAX, CABS (V2) )  .GT. 0. )  JCaORRECT-1 
M-INT(XMAX/  (CABS  (VI) ) )  +  NUCRRECT 
N2«INT(XMAX/  (CABS(V2) ) )  +  K2  CORRECT 
I2-INT  (MAX/  (CABS  (V2) ) ) 


-A 


,*  j 


TO  CHECK  OUT  WHICH  REGION  ARE  THE  VECTORS  IN 

VECTOR- Vina  SIGN  +  V2*N2SIGN 
IF  ( (REAL  (VECTOR)  *AIWG  (VECTOR) )  .CT.O)  SAW-. TRUE. 

IF  ( (REAL  (VECTOR)  *A1WG  (VECTOR) )  .LT.O)  SAW-. FALSE. 
IF  ((REAL(VECTOR)*AIWyG  (VECTOR))  .W.O)  00  TO  3 

vector-2.  nanusiGN  +  v2*n2sign 

IF  ((REAL  (VECTOR) ‘AIMAG  (VECTOR)). OT.  0)  SAW-.TKJE. 

IF  ( (REAL (VECTOR) ‘AIMAG (VECTOR) )  .LT.O)  SAW-. FALSE. 
3  CHECX  K). 

CO  5  j-i.m 

DO  10  I-I2-N2 

IF  (CHECK. GE.  2.)  00  TO  50 


•  • 
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VECTOR»J*V1*N1SIGN  +I*V2*N2SIGN 
MM3-CABS  (VECTOR) 

MAKE  SURE  THAT  (BECK  IS  INCREMENTED  PROPERLY 

IF  ((NOT  SAJE)  .Mffi.  ( (REAL (VECTOR)  *Aim3 (VECTOR) )  .GT.O) ) 
,  CHECK  <HECX  +  1. 

IF  ((SAME) .AND.  ( (REAL  (VECTOR)  *AIMMG  (VECTOR) )  .LT.O) ) 

.  CHECK -CHECH  +  1 

MAKE  SURE  THE  POINT  GENERATED  IS  INSIDE  THE  CIRCLE 

IF  (MM3.GT.Amx)  <30  TD  5 

MAKE  SURE  THAT  CHECK  IS  DECREMENTED  PROPERLY 

IF  ((NOT  SAMS ). AND.  ((REAL (VECTOR) *AimG (VECTOR) ).GT.O)) 
.  CHECK  KHECK  -  1. 

IF  ( (SAME)  .AND.  ( (REAL  (VECTOR)  *AimG  (VECTOR) )  .LT.O) ) 
i  CHECK  <HECX  -  1. 

PHASE-ATAN2  (AIMG  (VECTOR) ,  REAL  (VECTOR)  > 

CALL  NORMAL (PHASE) 

VECTDR-CMFLX  (PHASE. MM3) 

LOOSEARCH  (VECTOR,  DFILE.MSIZE.M) 

RJ-FLQAT(J) 

RI-FLQAT(I) 

IF  ( (HI SIGN.  LT.  0)  .AND.  (N2SIGN.  LT.  0) ) 

,  CALL  INSERT  (LOC.M.DFILE.  MS  IZE,  VECTOR. -RJ.-RI) 

IF  <  (NL SIGN. LT.O)  .AND.  (N2SIGN.GT.0) ) 

>  CALL  INSERT  (LDC.M.DfTLE,  MS  IZE,  VECTOR, -RJ.RI) 

IF  ( (NlSIGN.GT. 0)  .AND.  (N2SIGN.LT.0) ) 

.  CALL  INSERT  <IiX,M,DFILE,  MS  IZE,  VECTOR,  RJ.-RI) 

IF  ( (NlSIGN.GT. 0)  .AND.  (N2SIGN.GT.0) ) 
i  CALL  INSERT  (IOC, M.DTILE.MSIZE.VECIDR.RJ.RI) 

CONTINUE 
CONTINUE 
IlKJ-2 

DO  15  I-I2.N2 

DO  20  J-Il.Nl 

VECTORS  *V1*N1SIGN  +I*V2*N2SIGN 
MG-CABS  (VECTOR) 

IF  (MM3.GT.AmX)  GO  TO  15 

PHASE- AT AN2  (AOBG  (VECTOR)  .REAL  (VECTOR) ) 

CALL  NDREAL  (PHASE) 

VEdDR-CMPLX  (PHASE,  MM3) 

LOOSEARCH  (VECTOR, DFILE.MSIZE.M) 

RJ-FLQAT ( J) 

RI-FLQAT(I) 

IF  ( (NlSIGN.LT. 0)  .AND.  (N2SIGN.LT.0) ) 
i  CALL  INSERT  (LOC,M,DFILE,MSIZE,VECTOR.-RJ,-RI) 

IF  ((NLSIGN.LT.O)  .AND.  (N2SIGN.GT.0) ) 


I 


c 

C  FUNCTION:  SEARCH 

C 

C  THIS  FUNCTION  WILL  SEARCH  THE  LOCATION  WHERE  THE  ELEMENT 
C  FITS  INTO  A  DATA  ARRAY.  THE  ARRAY  ASSUFES  AT  LEAST  TOO 
C  MEMBERS.  THE  ELEMENT  MAY  HAVE  3  POSSIBILITIES :  - 
C  1)  THE  REAL  FART  IS  THE  SAFE  IN  THE  RETURNED  LOCATION 

C  2)  THE  REAL  PART  IS  THE  SAFE  IN  THE  RETURNED  LOCATION  +1 

C  3)  THE  REAL  FART  IS  BETWEEN  THE  ABCVE  TOO 

C 

C  THIS  FUNCTION  EMPLOYES  BINARY  SEARCH  TECHNIQUE  TD  LOCATE 
C 

C  ELEMENT  :  COMPLEX  ELEMENT  TD  BE  PLACED  WITH  PRIORITY  OF  THE 
C  REAL  PART  OVER  THE  DAGINARY  FART 

C  DFILE  :  ARRAY  FILE  13  BE  SEARCHED  AND  INSERTED 

C  MSIZE  :SIZE  OF  DFILE 

C  LAST  :NUFBER  OF  ELEMENT  IN  DFILE 
C 

FUNCTION  ffiARCH (ELEMENT, DFILE, MSIZE, LAST) 

COMPLEX  DFILE  (MSIZE,  2)  .ELEMENT 

LOC2-IAST 

LOC1-1 

IF  (REAL  (ELEMQiT) .  LT.  REAL  (DFILE  (IDC1 ,1) ) )  GO  TD  147 
IF  (REAL  (ELEMENT)  .GT.  REAL  (DFILE  (IOC2.1) ) )  GO  TD  140 
120  LOCM  -  (LOCI  +LOC2)  /2 

IF (LOCM. ED.  LOCI)  GO  TO  150 

IF  (REAL(ELEMENT).LT.REAL(DFILE(IOCM,l>))  GO  TO  125 
IF  (REAL  (ELEMENT)  .GT.  REAL  (DFILE  (LOCM.l ) ) )  GO  TO  130 
GO  TO  145 
125  IOC2  -LOCM 
GO  TO  120 
130  IOC1-UXM 
GO  TO  120 

140  SEARCH  -  L0C2 
GO  TO  155 

145  ffiARCH  -  LOCM 
GO  TO  155 

147  SEARCH  -  LOC1-1 
GO  TO  155 

150  SEARCH  -  LOCI 
155  RETUWJ 
am 
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SUBROUTINE: 


mawr 


THIS  SUBROUTINE  WILL  INSERT  THE  ELEMENT  INTO  THE  DFILE 

IF  THE  ELEMENT  HAS  A  DIFFERENT  PHASE  OR  A  SMALLER 

MAGNITUDE.  THE  PHASE  HAS  A  SENSITIVITY  SPECIFIED  BY  THE  ERROR. 

NSTART  :  LOCATION  OF  THE  ELEMENT 
LAST  .-NUMBER  OF  ELEMENTS  IN  DFILE 
DFILE  :  ARRAY  FILE  TO  BE  INSERTED 

MSIZE  :DI PENSION  OF  DFILE  _ 

ELEMENT  :  COMPLEX  ELEMENT  TO  BE  INSERTED 
VI INDEX  jNUPBER  OF  VI  USED 
V2  INDEX  : NUMBER  OF  V2  USED 

SUBROUTINE  INSERT  (NSTART,  LAST,  DFILE,  MSIZE,  ELEMENT,  VI  INDEX, V2 INDEX) 
COMPLEX  DFILE  (MSIZE.  2)  .ELEMENT 
PI-3.14159265 
ERROR-  0-01*PI/180. 

IF  (NSTART. NE. LAST)  GO  TO  200 

IF  (ABS  (REAL  (ELEMENT) -REAL  (DFILE  (PAST,  1) ) )  .LT.  ERROR)  GO  TO  290 

LAST  -  LAST  +1 

DFILE (LAST, 1 ) -ELEMENT 

DFILE (LAST, 2)  -CMPLX (Vl INDEX, V2 INDEX) 

GO  TO  295 

IF  (NSTART.  ED. 0)  GO  TO  202 

IF  (ABS  (REAL  (ELEMENT)  -REAL  (DFILE  (NSTART.l)  )  )  .  LT.  ERROR)  GO  TO  285 
IF  (ABS (REAL (DFILE (NSTARTf  1,1)) -REAL (ELEMENT) ).LT. ERROR)  GO  TO  280 
CALL  PUSH  (NSTARTH-1  .LAST, DFILE. MSIZE) 

DFILE  (NSTART+1 .1 )  -ELEMENT 

DFILE  (NSTARTH . 2 )  <MPLX  (VI  INDEX,  V2  INDEX) 

GO  TO  295 

IF  (AIMBG  (ELEMENT)  .GE.AIMVG  (DFILE  (NSTARTfl  ,1) ) )  GO  TO  295 

DFILE  (NSTARTU.l)  -  ELEMENT 

DFILE  (NSTARTU. 2)  KMPLX  (VI  INDEX,  V2  INDEX) 

GO  TO  295 

IF  (AH0G  (ELEMENT)  .GE.AIMAG  (DFILE  (NSTART,  1)))  GO  TO  295 

DFILE  (NSTART,  1 )  -ELEMENT 

DFILE  (NSTART,  2) -CMPLX  (VI  INDEX,  V2 INDEX) 

GO  TO  295 

IF  (AIMM3  (ELEMENT)  .GE.  AIMBG  (DFILE  (LAST,  1 ) ) )  GO  TO  295 

DFILE  (LAST.  1 )  -ELEMENT 

DFILE  (LAST,  2)  KMPLX  (VI  INDEX,  V2 INDEX) 

RETORN 
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SUBROUTINE: 


HJSH 


OHIS  SUBROUTINE  WILL  HJSH  AIL  THE  DATA  IN  'MTLE'  BY  1  POSITION 
ME  OPERATION  IS  SPECIFIED  BY  NSTART1  -TOE  BEGINNING  OF  HJSH 

LASTl  -TOE  LAST  LOCATION  OF  DETLE 

SUBROUTINE  HJSH (NSTARH .IAST1  ,DFTLE.MSIZE) 

COMPLEX  DFILE  (MSIZE,2)  ,HMPl ,HMP2 
JMASH-  NSTART1  +1 
DO  5  I-l.N 

TEMPI -DFILE (LASTl +1-1,1) 

TEHP2-DFILE  (LASTl+l-1, 2) 

DETLE  (LASTl +2-1,1 )  -TEMPI 
DETLE  (LASH +2-1,  2)  -TEMP2 
5  CONTINUE 

LASTl  -  LASH  +1 

RETUWJ 

END 


i 


I 


i 

I 


i  " 

B 


L 


c 

C  PROGRAM:  INTERPOL 

C 

C  THIS  EWE  RAM  WILL  INTERPOLATE  FOR  3  OPTIONS. 

C  FILE  FNAME  STORED  TOE  GRID  POINT  INFORMATION  GENERATED  EREVXDUSLY. 

C  THIS  EROGRAM  WILL  REPEAT  ITSELF,  EXIT  ONLY  BY  CONTROL  C 

C  THIS  FROG  RAM  LINKS  WITH  'PLOTLIB 

C 

C  THIS  PROGRAM  REQUIRES  AN  INRJT  OF  A  MEASUREMENT  FILE  :FNAfC 
C  THE  STRUCTURE  OF  ENAME: 

C  FMIN.FMAX  -MAXIMUM  AND  MINIMJM  FREQUENCY  USED  IN  HZ(2E15.8) 

C  EDELTA  -NORMALIZATION  FACTOR  IN  ERBQUENCY(1E15.8) 

C  NPOINT  -NUMBER  OF  POINTS  IN  THIS  FILE  (18) 

C  ISO, RADIUS- ISO' T'  IF  ISOTROPIC  SAMPLING  WAS  DONE(A2) 

C  -ISO'F'  IF  ISOTROPIC  SAMFLDC  WAS  NOT  DONE 

C  -RADIUS  IS  TOE  ISOTROPIC  CIRCLE  RADIUS  USED  IN  SECOND (E15. 8) 

C  V1.V2  -VECTORl  AND  VECIOR2  IN  FREQUENCY  D0MAIN(4E15.8) 

C  UI.U2  -VECTORl  AND  VECTDR2  IN  TUC/SPACE  EOMAIN(4E15-8) 

C  *.*.*, ‘-REAL  AND  IMAG  MASUREMENTjNl ,N2  FDR  TOE  *ASUREMENT(4E15.8) 

C  (MAKE  SURE  ALL  TOE  ABOVE  ARGURMENTS  ARE  SEPARATED  BY  COMMAS  IN  FNAME  1) 
C 

C  NL*V1,N2*V2  SPECIFIES  TOE  SAMPLING  LATTICE 

C 

C  LINK  INTERPOL . WFILE ,  SUM,  CDNST1 .  CDNST2 .  FUNC1 .  FUNC2 .  EUNC3 .  SINC,  - 
C  RHOT.ELIM,  'PLOTLIB 

C 

COMPLEX  SUM, DATA (360)  ,V1  .V2. U1.U2. FRED, CDEFFdOOOO. 2)  .TEMPI 
CHARACTER  ISO 
CHARACTER*10  ENA* 

MSIZE-10000 

C 

C  NDATA  sSIECTFIES  TOE  NUMBER  OF  ANGULAR  POINTS 

C  HIM  :SEECIfTES  TOE  NUEBER  OF  FREDUENCY  POINTS 

C 

NDATAf360 

NUfM=201 

PI-3.14159265 

WRITE (6.*)  '  FILE  WHICH  HAS  GRID  POINT  MEASUREMENT?' 

ACCEPT  50.  FNAtC 
50  FORMAT  (A10) 

OPEN  (UNIT-9.  NAME-FNAME, TYPE- 'OLD') 

C 

C  READ  IN  NECESSARY  DATA 

C 

READ(9,550)  FMIN.FMAX 

550  FORMAT (2E15.8) 

READ(9,551)  {DELTA 

551  FORMAT  (E115. 8) 

R£AD(9,552)  NPOINr 

552  FORMAT  (18) 

IF  (NFDINT.GT.MSIZE)  GO  TD  8888 


f 


153 


r 


READ(9,553)  ISO, RADIUS 
FORMAT  (A2- El  5. 8) 

READ(9.554)  V1.V2 
FORMAT  ( 2E1 5  -  8 , 2E1 5 . 8) 

R£AD<9.554)  U1.U2 
VI =2 . *PI *FDELTA*V1 
V2=2 .  *PI  *FDELTA*V2 
U1=UL/  (FDELTA) 

U2=U2/FDELTA 

OPTIONS  TD  GUARD  AGAINST  FEASUREHENT  DATA  IN  FORMAT 
OTHER  THAN  REAL  AND  LFftGINERY 

WRITE (6,*)  'DATA  IN  WHICH  EORMAT:l)  REAL  AND  IFftG’ 

WRITE (6.*)  '  2)  AMP (LINEAR)  AND  PHASE (RADIAN) 1 

WRITE  (6,*)  '  3)  AMP(DB)  AND  PHASE  (RADIAN)  ' 

WRITE  (6,*)  1  4)  AMP(DB)  AND  PHASE  (DEGREE)  ' 

ACCEPT  *.NTDM1 

IF  ((NCDM1.GT.4)  .OR.  (NCDMl.LT.l) )  GO  TD  570 
DO  560  I=l.NFOINT 

READ (9- 554)  XREAL.XIMAG,CDEFF(I,2) 

IF(NCDMl  .00.1)  03  ID  580 

AMP=XREAL 

FHASE^XIMRG 

IF((NCDM1.E0.3).OR.  (N03M1.ED.4))  AMP=10**  (AMP/20) 
IFtNCDMl-ED. 4)  PHASE=IHASE*PI/180. 

XREAL=AMP*CDS  (PHASE) 

XIMAG=AMP*S  IN  ( EHASE ) 

OOEFF  ( 1 , 1 )  =<J1HJC  (XREAL  -  XIMAG) 

CONTINUE 

CLOSE (UNIT=9 ,DISP= ' SAVE' ) 

ASK  FOR  THE  TYPE  OF  INTERPOLATION 

WRITE (6,*)  '  3  CHOICES  :1)  INTERPOLATE  FOR  1  ASPECT  ANGLE  AND  OUTPUT 
WRITE (6.*)  '  201  FREQUENCY  POINTS.  ' 

WRITE  (6.*)  '  2)  INTERPOLATE  FOR  1  FREQUENCY  AND  OUTPUT  ' 

WRITE (6.*)  '  360  DEGREE  JOINTS.’ 

WRITE (6.*)  '  3)  INTERPOLATE  INTO  2-0  SQUARE  GRID’ 

ACCEPT  *,  NQ3M 

THE  BRANCHING  OF  DECISION 

IF  (NCDM.BD.2)  GO  TD  200 

IF  (NCDM.EO.  3)  GO  TD  300 

IF  (NCDM.NE.1)  GO  TD  5 

TD  INTERPOLATE  FOR  1  ASPECT  ANGLE  BUT  201  FREQUENCY  POINTS 

WRITE (6.*)  'FREQUENCY  RANGE?  LOW  TD  HIGH  IN  HZ' 

ACCEPT  *.  FMINS,  FMAXS 


IF  (FMINS.GE.  FWAXS)  00  TD  10 

WRITE (6,*)  '  WHICH  ASPECT  ANGLE? (degree) ' 

ACCEPT  *.  ANGLE 
ANGLE=ANGLE*PI/180. 

FINCR= (FMAXS-FMINS) /(NJM-1) 

DO  15  1=1, NUM 

WRITE(6,*)  'WORKING  CN:  '.I,'  IAST  ONE:  '  .NUM 
APREQ-FMINS  +  FINCR*  (I — 1) 

IF  (A5REQ.EQ.0.)  00  TD  13 
C 

C  MAKE  SURE  THE  INTERPOLATION  IS  DONE  BETWEEN  THE  MEASURED  FREQUENCY 
C  RANGE 
C 

IF  (((AFREQ)  .GT.EMAX)  .OR.  ( (AFREQ)  .LT.  EMIN) )  00  TD  13 
FRED=CMPLX(AfREO*ODS  (ANGLE) , APRED*SIN  (ANGLE) ) 

DATA(I)  =  SUM  (FREQ,  CD  EFF.MSIZE,  VI  .V2.UL.U2.  ISO.  NIOINT,  RADIUS) 
00  TD  15 

13  DATA(I)-(0.,0.) 

15  CONTINUE 

C 

C  PLOT  AND  WRITE  FILE 

C 

CALL  RHCT  (DATA,  NJM.  FMINS,  FINCR) 

CALL  WFILE  (DATA.  NDATA,  NUM) 

QO  TD  5 
C 

C  SECOND  SECTION  OF  THE  PROGRAM 

C  TD  INTERPOLATE  ' FOR  1  FREQUENCY,  BUT  360  DEGREE  POINTS 

C 

200  WRITE (6.*)  ’WHICH  FREQUENCY? (HZ)  1 
ACCEPT  *,  AFREQ 
DO  205  1*1 .NDATA 

WRTTE(6.*)  'WORKING  CN:  '.I,'  LAST  ONE:  '.MATA 
ANGLE>(I-1)*PI/180. 

IF  (APRBQ.EQ .0.)  GO  TD  203 
C 

C  MAKE  SURE  THE  INTERPOLATION  IS  DONE  BETWEEN  THE  MEASURED  FRBOUENCY 
C  RANGE 

C 

IF  ( ( (AFREQ)  .GT.  EMAX)  .OR.  ( (AFREQ)  .LT.  FMIN) )  00  TD  203 
FRED-CMFLX  (AEREQ*CDS  (ANGLE)  ,AFRBQ*SIN (ANGLE) ) 

DATA  (I)-  SUM(FREQ,CDEFF,MSIZE,V1  ,V2  .U1  .U2.  ISO,  NFOD/T,  RADIUS) 
00  TD  205 

203  DATA(I)*(0.,0-) 

205  CONTINUE 

C 

C  PLOT  AND  WRITE  THE  FILE 

C  FOLARP  :A  SUBROUTINE  IN  '  PLOTLIB 

C 

CALL  FOLARP  (DATA.  3. 5. 3,1  .MATA.1,1) 

CALL  WFILE  (DATA.  NDATA.  NDATA) 


o  n  n  n  n  o 


00  TD  5 

THU®  PART  OF  THE  PROGRAM 
TD  DnnUOLAIE  INTO  A  2-D  SOUWE  GRID 

(DOE  TD  TIEE  CONSUMPTION.  THIS  IS  NOT  IMPLEMENTED.  BUT  IT  CAN 
BE  DONE  SIMILAR ILY  AS  FART  1  AND  2) 


300  WRTTE(6.*)  'NOT  READY  YETI' 

00  TD  5 

8888  WRITE (6.*)  'ERROR: THE  SIZE  OF  INPUT  FILE  IS  LARGER  THAN  SFECLFTED! 
00  TD  5 


COMPLEX  FUNCTIONS  SUM 

THIS  FUNCTION  WILL  DO  THE  SUMMATION  OF  F  (Nl*Vl+N2*V2)  *G  <W-ta*Vl-N2*V2) 
(BOUATIDN  2-9  IN  THIS  WRITE  UP) 

FOR  ALL  AVAILABLE  VALUES  OF  ML  AND  N2  IN  THE  FILE  FNAME 

FRBO  : COMPLEX  (F1.F2) ;  LOOiTIDN  IN  THE  FREQUENCY  PLANE 

CD  EFT  :  COMPLEX  ARRAY  STORING  THE  SAMPLED  VALUES 

MSI2E  .-SIZE  OF  CDEFF 

V1.V2  : VECTOR  1  AND  2  IN  THE  FREQUENCY  DOMAIN 

U1.U2  : VECTOR  1  AND  2  IN  THE  SPACE  DOMAIN 

ISO  :T»  ISOTROPIC  SAMPLING  ,  F-NDN-ISOTRQPIC  SAMPLING  (Al) 

NPOINT  sNUWffi  OF  NON- ZERO  ELEMENTS  IN  CDEFF 
RADIUS  :  RADIUS  OF  THE  ISOTROPIC  CELL 

THIS  REQUIRES  THE  SUPPORT  OF  CDNST1  .CDNST2 

COMPLEX  FUNCTION  SUM  (FRED. CDEFF.  MS  IZE.V1.V2.U1.U2.  ISO,  NPOINT,  RADIUS) 
CHARACTER  ISO 

COMPLEX  FREQ, VI . V2-U1.U2.W,  CDEFF  (MSIZE.  2) 

1*1*3.14159265 

SUM«CMPLX(0.,0.) 

CONVERSION  OF  FREQUENCY  INTO  RADIAN  EREQUENCY 

FRBC«.*PI*FREQ 
DO  5  I«1  .NPOINT 

W»  FREQ- REAL  (CD  EFT  (1, 2)  )*Vl-AIMWG(aOETF(I,2)  )*V2 
IF  (ISO.EQ. 'T')  SUM=SUM  +  CDEFF (1,1 ) *CDNST1  (RADIUS, W) 

IF  (ISO.EQ.  'F')  SUF^SUM  +  CDEFF  (1,1)  *OONST2(U1.U2.W) 

CONTINUE 

RETURN 


FUNCTION 


acNsn 


c 
c 

C  mis  FUNCTION  IS  FDR  RECONSTRUCTION  OF  ISOTROPIC  SAMPLING 
C  CALCULATING: 

C  (l/((R**2)«Wl*(Wl**2-3*W2**2)))*{  2*Wl*Q0S  (EWl/SQRT (3) )  *OOS (R*W2) 

C  -2*W1*CDS  (2*R*W1/SQRT<3) ) 

C  ~2*SQRT(3)*W2*SIN(R*W1/SQRT{3))*SIN(R*W2)} 

C 

C  mis  REQUIRES  THE  SUPPORT  OF  FUNQ.FUNC2.FUNC3.FLIM 

C 

C  R:  CONSTANT 

C  W:OOMILEX  (W1.W2) 

C 

FUNCTION  CONST!  (R.W) 

COMPLEX  W 
DIFENSIDN  ABO, 4) 

TOLER-1 - 
Wl-REAL(W) 

W2«AIMRG(W) 

OSORTO.) 

D1«(W1-C*W2)*R 
D2e  ( vfl.  +C*W2 )  *R 
C 

C  mE  RECONSTRUCTION  FUNCTION  WILL  BLOW  UP  WHEN  DVfl-0  OR 
C  2)WL**2=!3*W2**2 

C 

IF  ((ABS(Dl).GE.  TOLER).  AND.  (ABS(D2)  .GE.  TOLER))  GOTO  999 
IF  ( (Dl.ED.O.)  .OR.  (D2.BO.0.) )  00  TD  888 
C 

C  THERE  IS  CNE  (ORE  PROBLEM  DURING  IMPLEMENTATION: 

C  mE  COMPUTER'S  UNDERFLOW  PROBLEM. 

C  THEREFORE  FOR  TDLER-1.,  mE  FUNCTION  IS  INTERPOLATED 

C  WITH  THREE  POINT  MATCHING  POLYNOMIAL. 

C  mE  COEFFICIENTS  OF  mE  POLYNOMINAL  IS  CALCULATED  BY  GAUSSIAN 

C  ELIMINATION.  THIS  PROGRAM  IS  FURNISHED  BY  MR.  BING  KWAN 

C  ELIM(AB,3,4.3) 

C 

IF  ((ABS (02) ).LT. TOLER)  GO  TO  555 
AB(1.2)  -( (W1*R)+(TOLER)  )/C 
AB(2.2)»(W1*R/C) 

ABO. 2)  •( (WL*R)  -  (TOLER)  )/C 
GO  TD  777 

555  AB(1.2)  —  ( (TOLER) +(W1*R))/C 

AB(2-2)— W1*R/C 
ABO, 2)  *(  (TOLER)  -WL*R)/C 
777  AB(1.4)-FUNQ(R,W1.  (AB(1,2))/R) 

AB(2.4)-FUNC2(R,  (AB(2.2))/R) 

AB(3»4)«FUNQ(R,W1.  (AB(3,2))/R) 

DO  800  J-1,3 

AB(J,3)-1.0 


nnnnnnnnn 


FUNCTION 


0CNST2 


THIS  FUNCTION  IS  TD  RECONSTRUCT  RECTANSLUAR  SAMPLING 
CALCULATING: 

(SIN  ( .5*  (VI  .W) ) )  *  (SIN  ( .5*  CV2  .W) )  /  ( .25*  (VI  •W)  *  (V2  .W) ) 

THIS  FUNCTION  REQUIRES  THE  SUPPORT  OF  SINC 

FUNCTION  CCNST2 (VI .V2.W) 
complex  V1.V2-W 

Xl-REAL  (VI )  *REAL  (W)  +AIMG  (VI )  *AI>BG  ( W) 

X2-REAL  (V2)  *REAL  (W)  +AIMAG  (V2)  *AIMG  (W) 

C0NST2-SINC  (XI/ 2 . )  *S  INC  (X2/2 . ) 

RETURN 

END 
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SUBROUTINE 


RHOT 


THIS  SUBROUTINE  WILL  DO  RECTANGULAR  PLOT 

FOR  REAL  AND  IMAGINERY  HOT  OR  (OGNTIUDE  AND  PHASE  HOT 

DATA  : COMPLEX  ARRAY  TO  BE  PLOTTED  (SIZE-360) 

LAST  :LAST  ELEMENT  IN  THE  ARRAY 
FMIN  : SMALLEST  VALUE  ON  THE  AB SCISSOR 

FINCH  :  INCREMENT  ON  THE  AB  SCISSOR 


SUBROUTINE  RHOT (DATA, LAST, FMIN. FINCR) 

CHARACTER  (DM 
COMPLEX  DATA  (360) 

DIMENSION  YAXISl ( 201 ) , YAXIS2 ( 201 ) , XAXIS ( 201 ) 

WRITE (6,*)  'DO  YOU  WANT  REAL  AND  IMAGINERY  PL0T7Y/N' 
ACCEPT  6,  COM 
FORMAT  (Al) 

IF  (COM.  ED.  'N' )  GO  TO  15 
IF  (GOM.NE. '  Y' )  GO  TD  5 
DO  10  I-l.LAST 

YAXISl  (I )  -REAL  (DATA  (I ) ) 

YAXIS2 ( I ) =AIMAG  (DATA (I ) ) 

XAXIS(I)  »(I-1)  *FINCH  +  FMIN 
CONTINUE 
GO  TD  25 
DO  20  I-l.LAST 

XAXIS(I)«(I-1)*FIN®  +FMIN 
YAXISl ( I ) -CABS (DATA (I ) ) 

IF  (REAL  (DATA  (I)).  ED.  0.)  GO  TD  18 

YAXIS2 ( I ) -ATAN2  (AIMAG  (DATA (I ) ) , REAL (DATA (I))) 

GO  TD  20 

YAXIS2(IM). 

CONTINUE 

CALL  JLTPHG(XAXIS.YAXIS1,201.201.1.0.1) 

CALL  ELTPHG(XAXIS.YAXIS2.201. 201.1.0.1) 

RETUMJ 
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C  SUBROUTINE  :ELIM 

C 

C  THIS  SUBROUTINE  IS  COURTESY  OF  MR.  BING  KWANI16] 

C 

SUBROUTINE  ELIM CAB, N. NP, NDIM) 

DIMHSTON  AB(M>IM,NP) 

C  THIS  SUBROUTINE  SOLVES  A  SET  OF  LINEAR  EQUATIONS. 

C  THE  GAUSS  ELIMINATION  *THOD  IS  USED,  WITH  PARTIAL  PIVOTING. 

C  MULTIPLE  RIGHT  HAND  SIDES  ARE  PERMITTED,  THEY  SHOULD  BE  SUPPLIED 
C  AS  GOLUmS  THAT  AUG  MINT  THE  ©EFFICIENT  MATRIX. 

C  PARAMETERS  ARE: 

C  AB  COEFFICIENT  MATRIX  AUGMNTED  WITH  R.H.S.  VECTORS. 

C  N  ND.  OF  EQUATIONS 
C  NP  TOTAL  NO.  OF  COLUMNS  IN  AB 
C  NDIM  ND.  OF  ROWS  IN  AB 
C 

C  BEGINS  THE  REDUCTION 
NMl-N-1 
DO  35  i«i,wa 

C  FIND  THE  ROW  MJMBER  OF  THE  PIVOT  ROW.  WE  WILL  THEN 
C  INTERCHANGE  ROWS  TO  THE  PIVOT  ELEMENT  ON  THE  DIAGCNAL. 

IPVT-I 
IP1=I+1 
DO  10  J«IP1,N 

IF<ABS(AB(IPVT,I)).LT.ABS(AB(J,I)>)  IEVT-J 
10  CONTINUE 

C  CHECK  TO  BE  SURE  THE  PIVOT  ELEMENT  IS  NOT  TOO  SMALL,  IF  SO 
C  PRINT  A  MESSAGE  AND  RETURN. 

IF(ABS(AB(IIVT,I))  .LT.l.E-5)  GO  TO  99 
C  NOW  INTERCHANGE  ,  EXCEPT  IF  THE  PIVOT  ELEMENT  IS  ALREADY  ON 
C  THE  DIAGCNAL,  DO  NOT  NEED  TO. 

IF(IPVT.EQ.I)  GD  TO  25 
DO  20  JCDL“I,NP 
SAVE- AB(  I,  JCDL) 

AB  ( I ,  JCDL)  «AB  ( IEVT,  JCDL) 

ABQPVT,  JCDL) -SAVE 
20  CONTINUE 

C  NOW  REDUCE  ALL  ELEMENTS  BEIOW  THE  DIAGCNAL  IN  THE  I-TH  ROW. 

C  CHECK  FIRST  TO  SEE  IF  ZERO  ALREADY  PRESENT.  IF  SO 
C  CAN  SKIP  REDUCTION  EOR  THAT  ROW. 

25  DO  32  JH0W-IP1,N 

IF  (AB  (JROW,  I)  .EQ. 0.0)  GO  TO  32 
RATIO-AB  (JROW,  I) /AB  (I,  I) 

DO  30  KC0L-IP1,NP 

AB  (JROW,  ROOD  -AB  (JROW,  ROOD  -RATIO  *AB  (I,  ROOD 
30  CONTINUE 
32  CONTINUE 
35  CDNTIMJE 

C  WE  STILL  NEED  TO  CHECK  A(N,N)  FDR  SIZE. 

IF (ABS(AB(N,N)). LT.l.E-5)  GO  TO  99 


! 


i 


r 

v 

i  • 


165 


r. 


C  NOW  WE  BACK  SUBSTITUTE. 

NPl«Nfl 

IX)  50  KOOL=NPl,NP 
AB  (N,  ROOD  =AB  (N,  KDDL)  /AB  (N,  N) 

DO  45  J=2,N 
NVBL=NPl-J 
L-NVBL+1 

VALUE-AB (NVBL, KOOL) 

DO  40  K«L,N 

VALUED VALUE-AB ( NVBL, K)*AB(K, KOOL) 

40  CONTINUE 

AB (NVBL, KOOL) =VALUE/AB (NVBL, NVBL) 

45  CONTINUE 
50  CONTINUE 
RETURN 

C  MESSAGE  FOR  A  NEAR  SINGULAR  MATRIX. 

99  WRITE (66,100) 

100  FORMAT (/,' SOLUTION  NOT  FOSSIBE.  A  NEAR  0  PIVOT  ENCOUNTERED. ') 


n  o  n  n  n  nn  n  on 


PROGRAM  : 


INTEGFFT 


•mis  PROGRAM  WILL  READ  FROM  A  COMPLEX  DATA  FHE:ENAME1  (MSIZE) 

WHICH  IS  PART  OF  AN  INTERGRAND  FOR  INTEGRATION. 

THE  DATA  ARGUMENT  HAS  A  RANGE  OF  10.360)  . 

(I.E.  MSIZE  IS  USUALLY  A  MJLTIPLE  OF  360  <680) 

FOR  8192>MSIZE>512  CHANGE  M.NPl  AND  THE  ARRAY  SIZE  OF:  DATA1  .DATA2.S 
THE  OTHER  PART  OF  THE  INTEGRAND  IS  SPECIFED  IN 
FUNC1.  THE  INTEGRATION  IS  LIKE  A  CONVOLUTION. 

THEREFORE,  THE  INTEGRATION  IS  USING  FIT  APPROACH 

A  RECTANGULAR  TD  POLAR  COORDINATE  FILE  FUST  BE  SUPPLIED 
THIS  IS  FOR  THE  TIME  PLOT  WHICH  IS  FNAME3. 

THIS  FILE  IS  GENERATED  BY  PROGRAM :R£CTPOLAR 
FTJAME3  FORMAT: 

NUMBER  OF  POINT:  (*) 

RADIUS ,  ANGLE  (RADIAN)  ,  X-CDORDINATE ,  Y-CLORDINATE  ( 4E1 5-8) 

THE  ONLY  PART  THAT  IS  RELATED  TD  THE  MENSA'S  ACTUAL  INTEGRATION  IS 
SHOWN  AT  A  LATER  SET  OF  OOMFENT  BEFORE  TEMPI  IS  CALCULATED 

LINK  INUGFFT.APIDRI, EXPAND, PICK,  'SSP 

OOMFLEX  DATAl  (512)  ,DATA2(512)  ,FACIDR.K>LAR(10000)  ,RECT  (10000)  .DUMMY 
DIFENSIQN  S (128) 

CHARACTER*10  FNAME1 . FNAME2 . FNAME3 

CHARACTER  ONE, OLD 

M*=9 

NPl=512 

FOL-129 

NPOLAR=10000 

PI-3.14159265 

XNOR  IS  USED  TD  FORMALIZE  THE  AMP  OF  THE  OUTRJT 
XN0R-1E-9 

THIS  TTFE  FACTOR  IS  EACH  UNTT  OF  THE  TIFC  PLOT 

TIME-1/ (0-4921 2598E9*31 . ) 

WRITE (6,*)  'FILE  NAME  OF  MEASURED  DATA:' 

ACCEPT  9000.FNAME1 

WRITE (6.*)  'FREQUENCY  USED:  (HZ)  ' 

ACCEPT*.EREQ 

WRITE (6.*)  'NUFEER  OF  POINTS  IN  THE  FILE:' 

ACCEPT*. MSIZE 

IF  THE  FUFBER  OF  POINTS  IN  THE  FILE  IS  LESS  THAN  NR,, 

THEN  THE  DATA  MAY  BE  ONE  QUAERANT  DATA  ONLY 


IFCMSIZE.GT.NPl)  QO  ID  998 
IF  (MSIZE.GT.NH)  GO  ID  HO 
120  WRITE (6,*)  'IS  THIS  CUE  QUAERANT  DATA7Y/N' 

ACCEPT  9020, CUE 
IF  (ONE.EQ.  *N')  GO  ID  110 
IF  (ONE.NE.'Y')  GO  TD  120 

130  WRITE (6,*)  'CDMHEX  CONJUGATE  IN  2ND  AND  3RD  QUAERAOTS7Y/N' 
ACCEPT  9020, OLD 

IF  ((OED.NE.'Y').AND.  (OED.NE.  'N'))  GO  TD  130 
C 

C  TD  GUARD  AGAINST  FOSSIBLITY  OF  MEASUREMENT  IN  OTHER  FORMAT 
C 

110  WRITE (6,*)  'EORMAT  OF  DATA:  1)REAL  AND  IMAGINARY' 

WRITE  (6.*)  '  2)  AMP  (LINEAR)  AND  PHASE  (RADIAN)  ' 

WRITE  (6,*)  '  3)  AMP(LINEAR)  AND  PHASE  (DEGREE)  ' 

WRITE (6,*)  '  4) AMP(DB)  AND  MASE (RADIAN)  ' 

WRITE  (6.*)  '  5)AMP(DB)  AND  PHASE  (DEGREE) ' 

ACCEPT  *,NFORM 

IFC(NEORM.GT.5).OR.  (NIORM.LT.l))  GO  TD  110 
WRITE (6,*)  'STORAGE  FILE  NAME:' 

ACCEPT  9000.ENAME2 

WRITE (6.*)  'NAAE  OF  THE  RECTANGULAR  TD  POLAR  FILE: ' 

ACCEPT  9000.ENAME3 
C 

C  INITIALIZATION 

C 

DO  10  I-l.NPl 

DATAl <I)«(0.,0.) 

10  CONTINUE 

C 

C  READ  IN  MEASUREMENT  LATA 

C 

OIHN(UNrr»8,NAAE-=FNAMEl.TYPE>'OLD') 

DO  60  I-l.MSIZE 

READ (8,9030)  DUMMY 
CALL  APTDRI  (DUMMY,  NFDRM) 

DATAl  (I) -DUMMY 
60  CONTINUE 

OCSE  (UNIT-8  ,DISP"  ’  SAVE' ) 

C 

C  EXPAND  TD  FILL  THE  HALF  OF  THE  SPAN  OF  THE  FFT  REPETITIVE  UNIT 
C 

IF  (MSIZE.LE.  NFL)  CALL  EXPAND  (DATAl , DATA2 . NP1 , MS IZ E,  NFL) 

IF  (MSIZE.GT.NH,) CALL  EXPAND  (DATAl  ,rATA2,NPl.MSIZE,NPl) 

IF  (MSIZE.GT.NH)  QO  TD  170 
C 

C  CASE  OF  FULL  HANE  DATA  INPUT 

C 

IF  (ONE.  EO.  'N')  CALL  EXPAND  (DATAl,  DATA2.NPl,NH,NPl) 

IF  (ONE.EQ.  'N')  GO  TD  170 
C 
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C  FILL  IN  TOE  OTOIR  HALF  OF  TOE  SEAN  BY  COMILEX  CONJUGATE  OR 
C  THE  SANE 
C 

MSIZE=NFL 

DO  160  1=1  .MSIZE-2 

IF  (OEO. E0.  'N' )  DATA1  (I+MSIZE)  “CATAl  (MSIZE-I) 

IF  (ODD.  HQ.  'Y')  DATA1  (I+MSIZE)  =<JONJG(DATAl(MSIZE-I)> 

160  CONTINUE 

HSIZE=2*MSIZE-2 
DO  165  I“1,MSIZE 

IFtODD.EQ.  'N')  DATAl  (I+MSIZE)  “DATA1  (I) 

IF  (DDD.EQ. 'Y')  DATAl  ( I+MSIZE) =C0NJG(DATA1  (I) ) 

165  CONTINUE 

MSIZ  E=MSIZE*2 
IF(MSIZE.GT.NPl)  GO  TO  998 
170  MSIZE=NP1 
C 

C  DELTA  :SIZE  OF  EACH  ANGLE  INCREMENT 
C 

DELTA=2.*PI/MSIZE 

OPEN  (UNTP=8,NAfE=FNAME3  ,TYPE='OLD' ) 

READ  (8,*)  NFOINT 
IF  (NPOINT.CT.NEOLAR)  GO  TO  998 
C 

C  READ  TOE  RECTANGULAR  ID  POLAR  COORDINATE  CONVERSION  FILE 
C 

DO  100  I«l,NEOINT 

READ(8,9005)  FOLAR(I)  ,RECT(I) 

100  CONTINUE 

CLOSE  (UNnya.DISFta'SAVE’ ) 

CALL  EORTCDATAl.M.S.l.IFERR) 

C 

C  TOE  VALUE  OF  AMP1  —1  IS  CNLY  A  DUMMY  TO  START  TOE  ROUTINE 
C 

AMPl— 1. 

DO  200  J-l.NFOINT 

WRITE(6.*)  'WORKING  ON:'. J,'  LAST  ONE: '  .NPOINT 
IF  (REAL  (IOLNUJ)). HQ.  AMPl)  GO  TO  50 
AMPl “POLAR (J) 

DO  17  I“1,NP1 

OATA2(I)*(0.,0.) 

17  CONTINUE 

DO  20  1=1  .MSIZE 

TOETA“(I-1)*DELTA 

PHI“2.  *PI  *  (REAL  (POLAR  ( J) ) )  *TUE*FREQ*CDS  (TOETA) 
C 

C  FUNC1:  EXP(JWT) 

C  THEREFORE,  XREAL-CDS  (PHI) 

C  YIMAG“SIN(MI) 

C 


XREAL-CDS(PHI) 


YIWVG-SINtffil) 

CATA2(I)KXPLX(XR£AL,YI»G) 

20  CONTINUE 

C 

C  GO  TO  FREDUENCY  DOMAIN 

C 

CALL  TOFT(DATA2.M,S.2.IFERR> 

C 

c  convolution  in  tub  <->  multiplication  in  frequency 
c 

DO  25  I-l.NPl 

DATA2(I)=OATAl (I)*DATA2(I) 

25  CONTINUE 

C 

C  GO  BACK  TO  TOE  DOAMIN 

C 

CALL  EORT(DATA2.M,S.-2.IFE»R) 

C 

C  NOW  TOLAR  (J)  IS  THE  INTEGRAL  VALUE  OF  2D  FOURIER  TRANSFORM 
C  WITH  IMRJLSIVE  RADIUS  VALUE  (MENSA'S  INTEGRAL) 

C 

50  ANGLB-AIMAG  (TOLAR  ( J) ) 

TOLAR  ( J)  “PICK  (DATA2 .  NPl  .0 . ,  DELTA/  ANGLE) 

TOLAR ( J) “(POLAR ( J) /FLOAT (MSIZE) ) *FREQ*XNDR 

200  OONTMJE 
C 

C  WRITE  OUT  THE  FILE 

C  FILE  FORMAT: 

C  NTOINT:NUJBER  OF  JOINTS  (110) 

C  REAL , IWG INARY, X-CDORDINATE , Y-CDORDINATE  (4E15.8) 

C 

OPEN  (UNIT“9,NAW:“ENAM:2  .TYPE“ '  KEW' ) 

WRITE (9,9010)  NTOINT 
DO  300  I-1.NTOHTT 

WRITE  (9.9005)  POLAR(I)  ,RECT(I) 

300  CONTINUE 

CLOSE  (UNIT-9,DISP- '  SAVE' ) 

GO  TD  999 

998  WRITE(6,*)  'ERROR:  SIZE  OF  FILE  TOO  L/WGE1' 

C 

C  FORMAT  STATEMENTS 

C 

9000  FORMAT  (A10) 

9005  FORMAT (4E15. 8) 

9010  FORMAT  (110) 

9020  FORMAT  (Al) 

9030  FORMAT(2E15.8) 

C 


SUBROUTINE:  EXPAND 

THIS  SUBROUTINE  WILL  EXPAND  AN  ARRAY  DATA(MSIZE)  CONTAINING  'NFILL 
DATA  ELEMENTS  INTO  'NEXP'  DATA  ELEMENT  USING  LINEAR  INTERPOLATION 
2<NFUI<<NEXP 

WORK:  A  DUEMY  ARRAY  EOR  TEMPORARY  STORAGE 


SUBROUTINE  EXPAND  (DATA#  WORK*  MSIZE, NFILLf  NEXP) 
OOMH<  EX  DATA  (MSIZE), WORK  (MSIZE)  ,DIFF 
ERRDR*=lE!-6 

DELTAl-1 ./FLOAT (NF ELL-1) 

DELTA2-1 ./FLOAT (NEXP-1) 

DO  10  1-1. NFILL 

WORK(I)=OATA(I) 

OONTINUE 

DATA(1)=W0RK(1) 

DATA  (NEXP)  -WORK  (NFILL) 

DO  20  1*1 .NEXP-2 

X-(I)*DELTA2 
N“IOT  (X/DELTAl )  +1 
REMAEN=AMDD  Q( .  DELTA1  > 

RATIO* (X-DELTA1 ‘FLOAT (N-l) )  /DELTA1 
DIET*  (WORK  (ER-1)  -WORK  (N) )  *RATID 
DATA(I+l)=«ORK(N)  +DIFF 
CONTINUE 
RETURN 


COMPLEX  FUNCTION 


PICK 


C 
C 
C 

C  THIS  FUNCTION  WILL  RETURN  A  LINEAR  INTERPOLATED  COMPLEX 

C  VALUE  BACK. 

C  DATA:  ARRAY 

C  MSIZE:  SIZE  OF  THE  ARRAY 

C  FIRST:  FIRST/DELTA  IS  THE  FIRST  INDEX 
C  DELTA:  INCREMENT  OF  EACH  STEP  IN  THE  ARRAY 
C  VALUE:  VALUE/DELTA  IS  THE  LOCATION  IN  THE  ARRAY 
C  N  :  EXACT  LOCATION  OR  LOCATION  -1 

C  VARY  : AMOUNT  OF  ANGLE  DIFFERENCE 

C  FACTOR  : ADJUSTMENT  TD  DATA2(N) 

C 

COMPLEX  FUNCTION  PICK  (DATA* MSIZE,  FIRST, DELTA/ VALUE) 

COMPLEX  DATA  (MSIZE) /FACTOR 

IF  ( (VALUE.  LT.  FIRST)  .OR.  (VALUE.  GT.  (MSIZE*DELTA  +FTRST)))  GOTO  998 
N*INT{  (VALUE-FIRST)  /DELTA)  +1 
VARY- (VALUE) -( <N-1)*DELTA+FIRST) 

FACTOR-  ( (DATA(Nfl)  -DATA(N) )  /DELTA)  *VARY 
PICK-DATA  (N)  +FACTDR 
GO  TO  999 

998  WRITE  (6,*)  'ERHOR:  INTERPOLATED  POINT  IS  OUTSIDE  THE  RANGE  1 1 

999  RETURN 
END 
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PROGRAM  : 


RECTPOLAR 


ISIS  PROGRAM  WILL  CHANGE  RECTANGLAR  COORDINATES  ID 
POLAR  COORDINATES.  THE  PILE  'FNAME'  OUTHJTS  IN  THE 
FOLLOWING  FORMAT: 

RADIUS. IHI (RADIAN)  ,X.Y  (4E15.8) 

THE  FILE  WILL  ALSO  BE  ORGANIZED  FROM  SMALLEST  RADIUS 
TD  THE  LARGEST  RADIUS. 

NOTE  THE  MAXIMUM  X  COORDINATE  VALUE:  '(AX'  WILL  GIVE 

(MAX*2)**2  +2*2*MAX  +1  DATA  POINTS 

THEREFORE  A  GRID  OF  256X256  CAN  ACCOMODATE  WOC-127 

LINK  RECnOLAR,GEN.SEARCHl.INSERTl. PUSH. NORMAL 

COMPLEX  DATA (65536. 2) 

CHARACTER*10  FNAfC 

WRITE (6.*)  'FILE  NAtG  FOR  STORAGE  OF  RESULTS:1 
ACCEPT  5.FNAFE 
FORMAT (A10) 

WRITE (6,*)  'MUCIMUM  X  COORDINATE  VALUES:' 

ACCEPT  *,MAX 
MSIZE-65536 
PI-3.14159265 
(H) 

DATACM.IXKLXCO.  ,0.) 

DATACM,  2)  <MPLX(0.  ,0 .) 

THIS  DO  LOOP  GENERATES  ALL  THE  AXIS  POINTS 

DO  10  J-l.MAX 

Xl -FLOAT  (J) 

DATA  (M ,  1 )  <MPLX  (XI .  0 . ) 

DATA  (M.2XMPLX  Oil, 0.) 

I^fl 

DATA  (M,  1  >  -CMPLX  (XI .  PI  /2) 

DATACM.2)  -CHHJUO.  ,Xl) 

DATA  <M,  1 )  •OOLX  CXI .  PI ) 
nAIA(M,2)<MPLX(-X1.0.) 

j^fl 

DATA(M,lXMHJCQQ,3.*PI/2.) 

DATA(M,2)<MPLX<0.,-X1) 

CONTINUE 

WRITE (6.*)  'CALCULATING  FIRST  QUADRANT  FOIOTSl' 

CALL  GEN(MAX.DATA< MSIZE.M.l  .1) 

WRITE  (6,*)  'CALCULATING  SECOND  QUADRANT  10  INIS  l' 
CALL  GEN  (MAX.  DATA.  MS IZE.M. -1.1) 

WRITE (6.*)  'CALCULATING  THIRD  QUADRANT  POINTS  1' 


I 

►  «. 


I 


I 


E 

t  • 


* 


►V 

i 


& 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


10 

5 


SUBROUTINE:  GEN 

THIS  SUBROUTINE  WILL  GENERATE  OFF  AXIS  POINTS 
FOR  V1*(1,0) ,V2«(0.1) 

THIS  IS  A  MODIFIED  VERSION  OF  GEN1  IN  PIGRID 

MAX  :  MAXI  HIM  INDEX 
DFILE  :CX>MFLEX  ARRAY  FOR  STOAGE 
MSIZE  :DIHNSIDN  OF  DFILE 
M  :NUFBER  OF  ELEMENT 

NISIGN  :SKN  OF  Nl  (N1*V1) 

N2SIGN  :SIGN  OF  N2  (N2*V2) 

THIS  REQUIRES  THE  SUPPORT  OF  SEARCR1  .INSQtTl, NORMAL 

SUBROUTINE  Gm(MAX.DFTXE,MSIZE,M,Nl.SIGN.N2SIGN) 

EXTERNAL  SEAR  CHI 

OOMFLEX  DFILE  (MSIZE.  2)  .VECIOR.Vl  .V2 

REAL  MAG 

VKMPLX<1.,0.) 

V2«<MFLX(0.,1.) 

DO  5  1*1 .MAX 

WRITE (6.*)  I 
DO  10  J-l.MAX 

VECTDR-jna*NlSIGN  +I*V2*N2S2GN 

mat^apc 

PHASE* ATAN2  (AIMAG  (VECTOR)  .REAL  (VECTOR) ) 

CALL  tORFAL  (PHASE) 

VECTOR-CMPLX  (MAG.  PHASE) 

LOO  SEAR  CHI  ( VECTOR,  DFILE,  MSIZE.  M) 

RJ-FLQAT(J) 

RI-FLQAT(I) 

IF  ( (NISIGN.  LT.  0)  .AND.  (N2SIGN.  LT.  0) ) 

1  CALL  INSERTl  (IOC,M,DFTLE,MSIZE,VECTDR.-RJ,-RI) 

IF  ((NISIGN.LT.O)  .AND.  (tCSKN.GT.O) ) 

2  CALL  INSERTl  (LOC.M, DFILE, MSIZE. VBCTDR.-RJ.RI) 

IF  ((NISIGN.GT.O).AND.  (N2SKN.LT.0) ) 

3  CALL  INSERTl  (LOC.M, DFILE, MSIZE, VECTOR, RJ.-RI) 

IF  ((NISIGN.GT.O)  .AND.  (N2SIGN.GT.0) ) 

4  CALL  INSERTl (IOC,M,DFILE,MSIZE,VBCTDR,RJ,RI) 

CONTINUE 

CONTINUE 

RETUMJ 

END 
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FUNCTION 


SEARCH! 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


120 


125 

130 

140 

145 

147 

150 

155 


THIS  FUNCTION  WILL  SEARCH  THE  LOCATION  WHERE  THE  ELEMENT 
FITS  INTO  A  DATA  ARRAY.  THE  ARRAY  ASSUMES  AT  LEAST  TWO 
MEMBERS.  THE  ELEMENT  FSVY  HAVE  3  POSSIBILITIES 

1)  THE  REAL  BART  IS  GREATER  THAN  THE  RETURNED  LOCATION 

2)  THE  REAL  BART  IS  THE  SAFE 

BUT  THE  IMAGINARY  PART  >RETURNED  LOCATION  VALUES 

3)  THE  REAL  AND  IMAGINARY  FART  ARE  THE  SAFE 

AS  THE  RETURNED  LOCATION  VALUES 


ELEMENT  :<30MPLEX  ELEMENT  TO  BE  INSERTED 
DFILE  :  CDMH. EX  ARRAY  TO  BE  SEARCHED 

F1SIZE  :SIZE  OF  DFILE 

LAST  : NUMBER  OF  NON-ZERO  ELEMENT  IN  DFILE 


FUNCTION  SEARCH1  (ELEMENT, DFILE. MSIZE.  LAST) 

COMILEX  DFILE  (MSIZE, 2)  .ELEMENT 

LOG-LAST 

LOC1-1 

IF  (REAL (ELEMENT)  .LT.REAL (DFILE (L0C1,1)>)  GO  TD  147 
IF  (REAL (ELEMENT) .GT.REAL (DFILE (LOG ,1)))  GO  ID  140 
IF  ( (REAL (ELEMENT) .  ED. REAL (DFILE (LOCI ,1 ) ) ) .AND. 

1  (AIMAG(ELEMENT)  .LT.AIMAG(DFILE(L0C1,1) )) )  GO  TD  147 

IF  ( (REAL  (ELEMENT)  .EQ.  REAL  (DFILE  (IOG.l) ) )  .AND. 

1  (AJMAG  (ELEMENT)  .GT.AIMAG (DFILE (L0G,1) )) )  GO  TD  140 

LOCM  -  (Ida  +IOG)/2 
IF(IOCM.EO.IOCl)  GO  TO  150 

IF  (REAL  (ELEMENT).  LT.REAL  (DFILE  (IOCM,l)>)  GO  TD  125 

IF  (REAL (ELEMENT)  .GT.REAL (DFILE (IOQl,l)))  GO  TD  130 

IF  (AIMBG (ELEMENT). LT.AIMRG (DFILE (L0CM,1)>)  GO  TD  125 

IF  (AIMBG  (ELEMENT)  .GT.  AIMBG  (DFILE  (IOCM,  1 ) ) )  GO  TD  130 

GO  TD  145 

IOG  -IOCM 

GO  TD  120 

LOa-UXM 

GO  TD  120 

SEARCHl  -  LOG 

GO  TD  155 

SEARCHl  -  IOCM 

GO  TD  155 

SEARCHl  -  IOC1-1 

GO  TD  155 

SEARCHl  -  IOC1 

RETUHJ 


noon 


c 

C  SUBROUTINE:  INSERT! 

C 

C  THIS  SUBROUTINE  WILL  INSERT  THE  ELEMENT  INTO  THE  DFILE 
C  IF  THE  ELEMENT  HAS  A  DIFFERENT  REAL  AND/OR  IMAGINARY  PART. 

C  TOE  REAL  FART  HAS  A  SENSITIVITY  SEEQFIED  BY  THE  ERROR1. 

C  TOE  DAG  INARY  BART  HAS  A  SENSITIVITY  STOQFIED  BY  ERRCR2. 

C  TOE  HUDRITY  IS  REAL  FART  FIRST  THEN  DAG  INARY  PART. 

C 

C  NSTART  :  LOCATION  OF  TOE  ELEMENT 
C  LAST  :NUfBER  OF  ECN-ZEHO  ELEMENTS  IN  DFILE 
C  DFILE  : ARRAY  FILE  TD  BE  INSERTED 

C  MSIZE  :SIZE  OF  DFILE 

C  ELEMENT  :(DMTt,EX  ELEMENT  TD  BE  INSERTED 
C  VI  INDEX  sNUEBER  OF  VI  USED 

C  V2  INDEX  :NU(BER  OF  V2  USED 

C 

SUBROUTINE  INSERT!  (NSTART.  LAST,  DFILE,  MSIZE,  ELEMENT, VI  INDEX,  V2 INDEX) 
OOMELEX  DFILE  (MSIZE. 2)  .ELEMENT 
PI-3.14159265 
ERRORl-  l.E-6 
ERHOR2-O.01*PI/180. 

IF  (NSTART.  JC.IAST)  GO  TD  200 

CASE  OF  NSTART— LAST 

IF(ABS(REAL(ELEMENT)-REAL(DfTLE(LAST,l))).GT.ERRORl)  GO  TD  190 
IF(REAL(ELEMENT).GT.REAL(DFILE(LAST,1)))  GO  TD  190 
IF CABS(AIWVG (ELEMENT) -AIEBG (DFILE (LAST.l)))  .LT.ERROR2)  GD  TD  295 
190  LAST  -  LAST  +1 

DFILE (LAST, 1 ) -ELEMENT 

DFILE  (LAST,  2)  KMELX  (VI  INDEX,  V2INDEX) 

GO  TD  295 

200  IF  (NSTART.  EQ.O)  GO  TD  202 

C  IF  (ABS(REAL(ELEMENT) -REAL (DFILE (NSTART, l))).GT.ERRCfa)  GOTO  202 
IF  (REAL  (ELEMENT)  .GT.  REAL  (DFILE  (NSEART.l) ) )  GO  TD  202 
IF  (ABStAXMAG  (ELEMENT)  -AIE9C  (DFILE  (NSTART,  1 ) ) )  .LT.  ERROR2)  GOTO  295 
202  CALL  HJSH (NSTART+l .LAST, DFILE, MSIZE) 

DFILE  (NSTARTU  ,1 )  -ELEMENT 

DFILE  (NSTART+-1 .2)  KMPLX  (VI  INDEX,  V2  INDEX) 

295  RETURN 
END 
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PROGRAM: 


SUM3D 


THIS  PROGRAM  CAN  CO  THE  FOLLOWING: 

1)  READ  A  DATA  FILE 

2)  SUMMING  OF  DATA  FILES 

3)  PLOT  A  DATA  FILE  IN  A)  3-D  PICTURE  PLOT 

B)  CONTOUR  PLOT 

4)  2-D  FOURIER  TRANSFORM 

5)  MODIFIES  FREQUENCY  DATA  BY  A)  1/JW 

B) -1/W**2 

6)  COSINE  TAPERING  LOW  PASS  THE  FREQUENCY  DATA 

7)  HIGH  BASS  THE  FREQUENCY  DATA 

8)  WRITE  OUT  THE  FILE 

THE  INRJT  (OR  OUTEUT)  FILE  FORMAT: 

NPOINT:NUM3ER  OF  POINTS  IN  THE  FILE  (110) 

DATA,  X-CDORDINA3E,  Y-CDORDINATE  ( 4E1 5 . 8) 

THE  CONTOUR  PLOT  LINKS  WITH  NCAR  HOTTING  PACKAGES 
THE  LINKING  REQUIRED  IS: 

SCASSIGN  ERA2:  I NCAR2.  NCARHOT.  HOTL  IB]  NCARJLOTJLIB 
SASSIGN  ERA2: [NCAR2.NCARLIB]  NCAR_LIB 

SASSIGN  EEU12:  [NCAR2. NCARHOT]  NCAR_PLOT 

SASSIGN  DRA2:lNCAR2.NCARHiOT.CHROME]  NCAR__HOTLCHROME 
$ASSIGN  DRA2:  [NCAR2.  NCAR  HOT.  TEST]  NCAR_PLOT_TEST 

SASSIGN  CRA2-.tNCAR2.NCARPLOT.DOC]  NCWULOTLDOC 

SASSIGN  ERA2:  [NCAR2  .NCARHOT.  TRANSl  NCAR^HOTLTRANS 
SLINK  SSUM3D,DPLOT,TRAN2D, WEIGHT, PL0T3D,AITORI, 'SSP, ’ELOTLIB, 
ERA2:  [NCAR2. NCARHOT. HOTLIBlNCARCDNRE/LIB, NCARDASHC/LIB,- 
NCARGRAPP/LIB , NCARGRAFH/ LIB , [ NCAR2 .NCARLIB] UTTLITY/LIB 

CDMH.EX  DATA(3I,31) , DUMMY, APTORI,CT 

DIMENSION  ARRAY1 (31,31) ,ARRAY2 (31.31) ,ARRAYP (31,31) 

CHARACTER*10  FNAME1.FNAME2 

CHARAOTR  ADD. OMEGA,  SPACE,  LP,  HP 

m*=5 

MSIZE=2**N1-1 
PI»3. 14159265 
GJ=CMPLX(0.1) 

K>*1SIZ  E/2+1 

INITIALIZATION 

DO  5  1*1  .MSIZE 

DO  6  J-l, MSIZE 

DATA  (I,  J)  *■(().  ,0.) 

CONTINUE 

CONTINUE 

WRITE (6,*)  'DATA  FILE  NAFE:  ’ 

ACCEPT  9000,ENAME1 


WRITE (6,*)  'WiAT  IS  THE  DATA  EORMAT:  1) 

WRITE (6,*)  '  2)  AME 

WRITE (6,*)  *  3)  AME 

WRITE (6,*)  '  4)  AME 

WRITE  (6,*)  '  5)  AME 

ACCEPT  *,NCDM 

IF  { (NCDM.GT.5)  .OR.  (N0DM.LT.1) )  00  TO  110 


1)  REAL  AND  IMAG' 

2)  AMP(LINEAR)  AND  PHASE  (RADIAN) ' 

3)  AMP  (LINEAR)  AND  PHASE  (DECREE)  ' 

4)  AMP(DB)  AND  PHASE  (RADIAN)  ' 

5)  AMPffiB)  AND  PHASE  (DEGREE) 1 


READ  IN  THE  DATA 

OPEN  (UNIT-10 .  NAME-FNAFEl .  TYPE- '  OLD' ) 
READ(10,9010)  NPOINT 
DO  10  I=*l. NPOINT 

READ (10, 90 20)  DUfMY.Xl.Yl 
CALL  APIORI  (DUfM T,  NCDM) 
M=INT(X1)  -HC 
N-INT(Yl)  +K 

DATA(M,N)  “  DATA(M,N)  -HXIFMY 
OONTINUE 

CLOSE (UNIT-10 ,DISP“ ' SAVE' ) 


CALL  EPLOT (DATA,  MS IZ  E,  ARRAY1 , ARRAY2 . ARRAYP) 

ADD  ANOTHER  FILE 

WRITE (6,*)  '  ADD  ANOTHER  FILE?  Y/N' 

ACCEPT  9030,  ADD 
IF  (AH).  ED.  'Y')  00  TD  8 
IF  (AED.NE.'N')  00  TD  20 

WRITE  INTO  A  FILE 

WRITE (6,*)  '  WRITE  YOUR  DATA  INTO  A  FILE?  Y/N' 

ACCEPT  9030, KEEP 
IF  (KEEP.EQ.  *N' )  00  TD  30 
IF  (KEEP.NE.'Y')  00  TD  25 
WRITE  (6,*)  'FILE  NAIG: ' 

ACCEPT  9000.ENAME2 

OPEN  (UNrT-8,NAME=FNAME2.TYPE>',NEW,> 

WRITE (8,9010)  MSIZE**2 
DO  27  I-l.MSIZE 

WRITE  (8,9020) ,  (DATA(I,J)  ,  FLOAT  (I-K)  ,ELOAT(J-K)  ,J-1.MSIZE) 
OONTINUE 

CLOSE (UNIT-8 ,DISP«' SAVE') 


FOURIER  TRANSFORM 

WRITE (6,*)  '  00  TD  FREDUINCY  DOMAIN?  Y/N' 
ACCEPT  9030, OMEGA 


81 


IF  (OMEGA.  ED. 'N')  GO  ID  50 
IF  (OMEGA. NE.  'Y')  GO  ID  30 
CALL  1RAN2D(DATA,MSIZE,-1) 

CALL  DPLOT (DATA,  MSIZE,  ARRAY1.ARRAY2.ARRAYP) 

GO  ID  25 

INVERSE  FOURIER  1RANSEORM 

WRITE (6.*)  '  GO  ID  TIME  DOMAIN?  Y/N' 

ACCEPT  9030, TIME 
IF  (TIME.BQ. *N')  GO  ID  60 
IF  (TIME.NE.  'Y')  GO  ID  50 

MODIFY  ID  PREPARE  FOR  STEP  OR  RAMP  RESPONSE 

WRITE (6,*)  '  WANT  ID  MODIFY  DATA?  BY  1)1' 

WRITE  (6,*)  '  2)  1/JW' 

WRITE (6,*)  '  3) -1/W**2' 

ACCEPT  *.M3D 

IF  ( (MOD.  GT.  3 )  .OR.  (MOD.  LT.  1 ) )  GO  ID  120 
IF  (MOD.  ED.  1)  GO  ID  95 

IF  (MOD.  ED. 3)  GO  ID  130 

MODIFY  BY  1/JW 

DO  125  1=1  ,MSIZE 

DO  127  J=1.MSIZE 

W=2. *PI*SORT(FLQAT((I-MSIZE/2-l)**2+(J-MSIZ E/2-1)  **2) ) 
IFCW.NE.O.)  DATA(I,J)=DATA(I,J)/(CJ*W) 

IF  (W.ED.OJ  DATA(I,  J)  =(0.  ,0.) 

OONTOUE 
CONTINUE 
GO  ID  145 

MODIFY  BY  -1/W**2 

DO  135  I=1.MSIZE 

DO  140  J=1  .MSIZE 

W=2.*PI*S0FT  (FLOAT ( (I-MSIZ E/2-1) **2+(J-MSIZE/2-l) **2) ) 
IF  (W.NE.O.)DATA(I, J) »-DATA(I»  J) / (W**2) 

IF(W.ED.O.)  DATA(I, J)  =(0.,0.) 

CONTINUE 

CONTINUE 

CALL  DPL0T(DATA,MSIZE,ARRAY1,ARRAY2.ARRAYP) 

COSINE  TAPERING  LOW  PASS 

THIS  IS  BY  SPECIFYING  THE  NUM3ER  OF  HARMONICS 

WRITE (6,*)  'WANT  ID  DO  COSINE  L0W-EASS7Y/N' 

ACCEPT  9030, LP 

IF  (LP,  ED.  'N')  GO  ID  100 


AD-A1C2  593  SPACE -FREQUENCY  SAMPLING  CRITERIA  FOR  ELECTRONAQNETIC 

SCATTERING  OF  A  FINITE  OBJECT<U>  OHIO  STATE  UNIV 
COLUMBUS  ELECTROSCIENCE  LAB  F  V  FOK  AUO  85 
UNCLASSIFIED  ESL-71419G-11  NO0814-82-K-SB37  F/G  28714 


IF  (LP.JE.  'Y')  00  TO  95 

WRITE  (6,*)  'COT-OFF  ELEMENT  NUfBER: 1 

ACCEPT  MJCOT 

DO  200  I-l.MSIZE 

DO  300  J-l.MSIZE 

DATA  (I,  J)  -DATA  (I,  J)  "WEIGHT  (NCUT,  I,  J.MSIZE) 

300  CONTINUE 

200  CONTINUE 

CALL  CPLOTCDATA.MSIZE,  ARRAY1.ARRAY2.ARRAYP) 

C 

C  HIGH  PASS  FILTERING 

C  AGAIN  BY  sreamtc  THE  NUfBER  OF  HARMONICS 
C 

100  WRITE<6,*)  'HIGH  PASS  FILTERING 7Y/N1 
ACCEPT  9030. HP 
IF(HP.EQ.  'N')  00  TO  400 
IF(HP.tC.'Y')  00  TO  100 
WRITE  (6.*)  'CUT-OFF  ELEMENT  NUJBERs ' 

ACCEPT  *.NCUT 
DO  500  I-l.MSIZE 

DO  600  J-l.MSIZE 

RFRE1>30RT  <  (FLOAT  (I-MSIZ  E/2-1) )  **2  +  (FLOAT  (J-MSIZ  E/2-1 )  )**2) 
IF  (RFREO. LE. FLOAT (NCUT) )  DATA(I, J) -<0.,0.) 

600  CONTINUE 
500  OONTOUE 

CALL  EPLOT(DATA,KSIZE,ARRAY1.ARRAY2.ARRAYP) 

C 

C  DO  FOURIER  TRANSFORM  NOW 

C 

400  CALL  TRAN2D  (DATA,  MSIZ  E.  1 ) 

WRITE  (6.*)  'THE  DATA  ARE  IN  TUC  DOMAIN  NOW  I ' 

CALL  DPL0T(DATA,MSIZE.ARRAY1.ARRAY2.ARRAYP) 

00  TO  25 

9000  FORMAT  (AlO) 

9010  FDRMAT(IIO) 

9020  FORMATME15.8) 

9030  FORMAT (Al) 

60  STOP 

END 
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SUBROUTINE 


DPLOT 


•THIS  SUBROUTINE  WILL  PLOT  3-D  FDR  A  COMPLEX  ARRAY  DATA 
EITHER  IN  A  3-D  PICTURE  OR  A  CONTOUR  PLOT 
THE  CONTOUR  PLOT  LINKS  WITH  NCAR  PLOTTING  PACKAGE 

DATA  :2-D  COMPLEX  ARRAY  TO  BE  PLOTTED 
MSIZE  :SIZE  OF  DATA 

ARRAYl.ARRAY2.ARRAYP:DUMff  ARRAYS  WITH  SAME  DIMHSTONS  AS  DATA 

SUBROUTINE  DPLOT (DATA/ MSIZE, ARRAY1  > ARRAY2. ARRAYP) 

COMPLEX  DATA  (MSIZE,  MSIZE) 

DlfENSION  ARRAY1  (MSIZE,  MSIZE)  ,ARRAY2  (MSIZE,  MSIZE) ,  ARRAYP  (MSIZE,  MSIZE) 
CHARACTER  ACT, PLO, COM, COM2, 00M3,CDM4.0DM5. LINE, PI CX 
WRITE  (6.*)  'DO  YOU  WANT  TD  PL0T7Y/N' 

ACCEPT  9030, COM 
IF  (COM.BD.  'N')  GO  TD  9999 
IF  (OOM.NE.  *Y')  00  TD  10 

NORMALIZE  THE  WHOLE  DATA  FILE  BY  SOM)  FACTOR 

0  WRITE (6.*)  'DO  YOU  WANT  ID  NORMALIZE  WITH  A  FACT0R7Y/N' 

ACCEPT  9030, OOM5 
IF(OOMS.EO. 'N')  00  TD  510 
IF(OOM5.NE.'Y')  00  TD  500 
WRITE (6,*)  'FACTOR:' 

ACCEPT  *,  FACTOR 
GO  TD  20 
0  FACTDR“1. 

WRITE  (6.*)  '  PLOT  REAL  AND  IMAGINARY?  Y/N' 

ACCEPT  9030,  ACT 
IF  (ACT.  ED.  'N')  00  TD  50 
IF  (ACT.NE.'Y')  00  TD  20 
DO  40  I«l, MSIZE 

DO  30  J-l  .MSIZE 

ARRAY1  (I,  J)  -FACTOR*REAL  (DATA  (I ,  J)  > 

ARRAY2(I,  J)“FACTOR*AMG(DATA(I,  J) ) 

OONTINUE 
CONTINUE 
GO  TD  100 

CONVERSION  INTO  LINEAR  AMPLITUDE  AND  PHASE  (RADIAN) 

DO  70  1-1  .MSIZE 

DO  €0  J“l, MSIZE 

ARRAY1  (I ,  J)  -FACTDR*CABS  (DATA  (I ,  J) ) 

IF  ((REAL (DATA (I, J) )  .BQ.O.)  .AND. 

(ADAS (DATA (I, J)). ED. 0.))  00  TD  55 
ARRAY2(I,  J)  -ATAN2  (REAL (DATA ( I,  J) )  ,AIMG(DATA(I,  J) ) ) 

00  TD  60 


55  ARRAY2(I,J)*). 

60  OCNTINUE 

70  CENTINUE 

c 

c***  *********** 

c 

100  WRITE(6.*)  *1) 3D  PICTURE, 2)CENTDUR  ROT' 

ACCEPT  9030, PICK 

IF  ( (PICX.fE.  '2')  .AM).  (PXOC.1C.  '1 '))  00  X)  100 

101  WRITE (6,*)  'SLOT  AMFLITUDE/REALTY/N' 

ACCEPT  9030, COM4 

IF((DM4.B0.  'N')  00  ID  150 
IF(OOM4.fE.  'Y')  00  ID  101 
IF  (PICK.  ED.  '1')  GO  X)  102 
C 

C  CONTOUR  HOT 

C 

IF  (PICK.HO.  '2') CALL  E2CfJTRCARRAYl.HSIZE.MSnE) 

00  X)  150 
C 

C  PLOT  3  D  PICTURE 

C 

102  CALL  VR/OTS (0.0.0) 

120  WRITE (6.*)  ’HOT  X-LINE, Y-LINE, OR  B07H7X,Y,B' 

ACCEPT  9030, LINE 

IF  ((LINE.tC.  ’X')  .AND.  (LINE.  1C.  *Y')  .AM).  (LINE.NE.  'B') )  QC 
CALL  SETLLINES  3D  (LINE) 

CALL  SETLR0TA3TDN  3D(0) 

XMAX-ARRAYl(l.l) 

XMB«WRAY1<1,1> 

DO  300  I-l.NSIZE 

CO  310  J-l.MSIZE 

XMAX-AMAX1  (ARRAYl  (I,J)  ,XMAX) 

XMP^AKINl  (ARRAYl  <I,J)  ,XHN) 

ARRAYP (I, J) -ARRAYl (I, J) 

310  CONTINUE 

300  CONTINUE 

WRITE  (6.*)  'SOLE  OF  AMPLITUDC/REAL  PLOT: ' 

WRITE  (6,*)  '  (MAX*'  ,XMAX, '  KH*>',XMXN 
WRITE(6,*)  '  REGDMCM)  SOALEi '  .0.5/(XHAX-XKN) , ')  * 

ACCEPT  ‘.SCALE 

WRITE (6,*)  'WAT  IS  THE  VIEW  ANCLE  W.R.T.  X-Y  PLANE?' 
WRITE (6,*)  '  (RECDMENDING : 45) ' 

ACCEPT  *.VTEW 

CALL  HjOT_3IL5URFACE (ARRAYP, MSIZE, MS IZE.VTEK.O. ,0. .SCALE) 
<001  WRITE (6,*)  'IS  THIS  THE  LAST  ILOT7Y/N' 

ACCEPT  9030,  PLO 

IF  ((ILO.fC.  'Y')  .AID.  (ILO.fC.  'N'))  00  X)  4001 
IF  (PLO.  BO. 'N')  NSIGN-1 
IF  (FLO.  BO. 'Y')  NSTGfM 
CALL  FLCTtO.  ,0.,NSIGN*999) 


HOT  IMAGINARY  OR  ERASE 


C 
C 

150  WRTTE<6.*)  'HOT  THE  IRASE/IJWSINARY7Y/N' 

ACCEPT  9030, COM3 
IF  (COM3. ED.  ’N')  CD  ID  5000 
IF  (OOM3.1E.  'Y' )  00  ID  ISO 
IF  (PICK. BO. '1')  00  ID  157 
C 

C  CONTOUR  HOT 

C 

CALL  EZCNTR(ARRAY2.KSIZE,MSIZE) 

GD  ID  5000 
C 

C  3  0  PICTURE  HOT 

C 

157  CALL  VHOTS  (0.0.0) 

220  WRITE<6.*)  'HOT  X-LINE,  Y-LINE,  OR  BOW?X,Y,B' 

ACCEPT  9030. LINE 

IF  ((LINE.PG.  'X')  .AND.  (UNE.NE.  'Y')  .AND.  (LINE.NE. ' B* )  >  00  ID  220 
CALL  SET-LINES  3D  (LINE) 

CALL  SET-ROTATION  3D(0) 

C 

C  FIND  IRE  MINIUM  AND  MAXINJM  IN  THE  HOTTING  SET 
C 

XMUMWRAY2U.1) 

XMDMRRAY2U.1) 

DO  400  I-l.KSIZE 

DO  410  J-l.MSIZE 

XMAX*AMAXl  (ARRAY2(I,  J)  ,XNAX> 

XMIN-AMIN1  (ARRAY2(I,  J)  ,XMIN) 

ARRAYP  (I ,  J)  MJWRAY2  ( I ,  J) 

410  OONTUUE 

400  CONTINUE 

WRITE  (6,*)  'SCALE  OF  PHASE/IMG  INARY  PICT:' 

WRITE (6.*)  '  (MAX-*  .XMAX, '  MH*'.XMIN 
WRITE(6,*)  '  REOOMCND  SOLE* '  ,0.5/ OWAX-XWN) ,’) ' 

ACCEPT  ‘.SCALE 

WRITE (6,*)  'WHAT  IS  VIE  VIEW  ANGLE  W.R.T.  X-Y  PLANE?' 

WRITE (6,*)  '  (REQ0M£NDING:45) ' 

ACCEPT  ‘.VIEW 

CALL  HOT-3D  SURFACE  (ARRAYP,HSIZE,NSIZE,VIEW.0.,0., SCALE) 

4000  WRTTE(6.*)  'IS  THIS  WE  LAST  HOT7Y/N' 

ACCEPT  9030,  FLO 

IF  ((ILO.NB.  'Y')  .MO.  (HD.fC.  'N'))  GO  ®  4000 
IF  (PLO.BD.'N')  NSIGf^-1 
IF  (PLO.BQ.'Y')  NSIGf^l 
CALL  HOT(0.,0.,NSIGN*999) 

5000  WRITE  (6.*)  'HOT  AGAIN7Y/N' 

ACCEPT  9030.  COM2 


SUBROUTINE 


TRAN2D 


1 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 
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C 
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THIS  SUBRCUNITNE  MILL  KXJRIER  TRANSFORM  ON  A  COMPLEX  ARRAY  IN  2D  PLANE. 
THE  PLANE  IS  ASSUMED  TO  HAVE  THE  ORIGIN  AT  THE  CENTRE  OF  THE  PLOT. 

THE  PLANE  HAS  EQUAL  EOBER  OF  UNITS  CN  EACH  SIDE  OF  THE  THE  AXES. 

I.E.  THE  SIZE  OF  THE  ARRAY  (MSIZE)  IS  AN  ODD  NUICER. 

WORKl  IS  THE  WORKING  ARRAY  WITH  SIZE  2**N  TO  JUST  FIT  THE  DATA  ARRAY. 


DATA  :  COMPLEX  2-0  ARRAY  TD  BE  TRANSFORMED,  AND  RESULT  STORED 

MSIZE  :SIZE  OF  DATA 
IFSET  s+1.00  TD  TIME 

s-l.QD  TD  FREQUENCY 

SUBROUTINE  TRAN2D (DATA, MSIZE,  IFSET) 

COMPLEX  WORKl (32.32)  .DATA (MSIZE, MSIZE) 

INTEGER  INV(8)  ,m(3) 

DlfENSIDN  S(8) 

N»(MSIZE+l)/2 
MS  IZ  El-32 

mm(1)*5 
fft(2)  -5 
MM  (3) «0 

PLACE  THE  CENTRED  PICTURE  INTO  THE  FOUR  CORNERS 

DO  10  I-l.MSIZEl 

DO  5  J-l.MSIZEl 

WORKl  (I,J)-(0.,0.) 

CONTINUE 
CONTINUE 
DO  20  I-l.N 

DO  15  J*1.N 

WORKl  (I,  J)  -OATA(W-I-l.NfJ-l) 

CONTINUE 
CONTINUE 
DO  30  I-l.N-l 

DO  25  J-l.N-1 

WORKl  (MSIZEl-JHT+1  .MSIZE1-NKJ+1)  "DATA  (I,  J) 

CONTINUE 
CONTINUE 
DO  40  1-1, N 

DO  35  J-l.N-1 

WORKl  (MS IZ  El -W J+l ,  I )  -DATA  ( J,  I+N-l ) 

WORKl  (I,  J+MSIZE1-N+1)  -DATA (I+N-l  ,J) 

CONTINUE 

CONTINUE 

FOURIER  TRANSFORM 

CALL  HARM  (WCRK1, ML  INV.S.  IFSET,  IFERR) 


v. 


J 


c 

c 


KJT  THE  FOUR  CDFfrOS  BACK  ID  THE  OUTRE 


DO  120  I-l.N 

DD  115  J-1,N 

t*TA(N+I-l,N*J-l)  **CRXl  (1,  J) 

115  CONTINUE 

120  CONTINUE 

DO  130  I-l.N-1 

DO  125  J-l.N-l 

DATA  (I,  J)  HCRKl  (MSIZEl-W-I+1  ,MSIZE1-N*J+1) 
125  CONTINUE 

130  CONTINUE 

DO  140  I-1,N 

DO  135  J-l.N-1 

DATA(J,  I-t-N-DHCRKl  (MSIZEl-NfJ+1 .1) 

CATA(I+N-l,J)^CRKl(I,J+4BIZEl-»fl) 

135  CONTINUE 

140  CONTINUE 
RETUR1 
END 


•ft.  jJj  *  A  ♦ 


FUNCTION 


WEIGHT 


THIS  FUNCTION  WILL  CALCULATED  THE  COSINE  LOW  BASS  WEIGHTING 
WITH  THE  ASSUMPTION  THAT  THE  FREQUENCY  RESPONSE  DATA  IS  LOCATED 
IN  THE  MIDDLE  OF  THE  HOT 

NCUT-  CUT-OFF  NUFBHt  FOR  THE  LOW  PASS  FROM  THE  ONUR 
NX  -  X -COORDINATE  IN  THE  MSIZE  X  MSIZE 
NY  -  Y-COORDINATE  IN  THE  MSIZE  X  MSIZE 
MSIZE)-  SIZE  OF  THE  ARRAY 

FUNCTION  WEIGHT  (NCUT*  NXrNYr  MSIZE) 

CHARACTER  LP 

PI-3.14159265 

X  >*BS  (FLOAT  (NX-MSIZE/ 2-1) ) 

Y-ABS  (FLOAT  (NY-MSIZE/2-1) ) 

RADIUS-SORT  OC**2+Y**2) 

IF  (RADIUS.  LE.  FLOAT  (NCUT) )  WEIGHX) . 5*  (1 .  +ODS  (PI *RADIUS/FLQAT  (NCUT) ) ) 
IF  (RADIUS.  GT.  FLOAT  (NOJT) )  WEIGHT>0. 

RETURN 
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